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ABSTRACT 

The Universe is opaque to extragalactic very high-energy gamma rays (VHEGRs, E > lOOGeV) because they 
annihilate and pair produce on the extragalactic background light. The resulting ultra-relativistic pairs are com- 
monly assumed to lose energy primarily through inverse Compton scattering of cosmic microwave background 
photons, reprocessing the original emission from TeV to GeV energies. In Broderick et al. (2012, Paper I of 
this three paper series), we argued that this is not the case; powerful plasma instabilities driven by the highly 
anisotropic nature of the ultra-relativistic pair distribution provide a plausible way to dissipate the kinetic energy 
of the TeV-generated pairs locally, heating the intergalactic medium (IGM). Here, we explore the effect of this 
heating upon the thermal history of the IGM. We collate the observed extragalactic VHEGR sources to determine 
a local VHEGR heating rate. Given the pointed nature of VHEGR observations, we estimate the correction for the 
various selection effects using Fermi observations of high and intermediate peaked BL Lacs. As the extragalactic 
component of the local VHEGR flux is dominated by TeV blazars, we then estimate the evolution of the TeV 
blazar luminosity density by tying it to the well-observed quasar luminosity density, and producing a VHEGR 
heating rate as a function of redshift. This heating is relatively homogeneous for z < 4, but there is greater spa- 
tial variation at higher redshift (order unity at z '--^ 6) because of the reduced number of blazars that contribute to 
local heating. We show that this new heating process dominates photoheating in the low-redshift evolution of the 
IGM and calculate the effect of this heating in a one-zone model. As a consequence, the inclusion of TeV blazar 
heating qualitatively and quantitatively changes the structure and history of the IGM. Due to the homogeneous 
nature of the extragalactic background light, TeV blazars produce a uniform volumetric heating rate. This heating 
is sufficient to increase the temperature of the mean density IGM by nearly an order of magnitude, and at low 
densities by substantially more. It also naturally produces the inverted temperature-density relation inferred by 
recent observations of the high-redshift Lya forest, a feature that is difficult to reconcile with standard reionization 
models. Finally, we close with a discussion on the possibility of detecting this hot low-density IGM suggested by 
our model either directly or indirectly via the local Lya forest, the Comptonized cosmic microwave background, 
or free- free emission, but find that such measurements are currently not feasible. 

Subject headings: intergalactic medium - BL Lacertae objects: general - gamma rays: general - cosmology: 
theory - large-scale structure of Universe 



1. INTRODUCTION 

The Fermi satellite and ground based imaging atmo- 
spheric Cerenkov telescopes such as H.E.S.S., MAGIC, 
and VERJTA^] have demonstrated that the ultra-high en- 
ergy Universe is teeming with energetic very high-energy 
gamma-ray (VHEGR, E > 100 GeV) sources, the extragalac- 
tic component of which mainly consists of TeV blazars 
with a minority population of other sources such as ra- 
dio galaxies and starburst galaxies. These VHEGR obser- 
vations are being used to constrain t he sites a nd mecha 
nisms of particle acceleration (|see, e. g .. iPaglione 

Domingo-Santamarfa & Torres 2005; Thompson ^ 

Persic et al. 2008; de Cea del Pozo et al. 2009; Repha eU et all 
2010t iLacki etal.l I201CI|I). dynamics of black hole jets 
[Tl974t IC ■ 



; et al.l 
I et al.' 



1996 



2007 



(see, e.g.. 



1989 



iJonesetal. 119741; iGhisellini&Marasch 
Ghisellini & Tavecchi _ol l2008t ll^vecchio & GhiseUini 
GhiselUni et al .1 120091) . and intergalactic magnetic fields 



2008 



' High Energy Stereoscopic System, Major Atmospheric Gamma Imaging 
Cerenkov Telescope, Very Energetic Radiation Imaging Telescope Amy Sys- 



riGMF: iNeronov & VovkllMot iTavecchio et al."201() 
Dermeretal.1 1201 U [Taylor et all 120111: iDolag et all 



2011 



2011 



Takahashi et al.ll2012l:1\Wk et al.ll2012l) 

While these objects have an interesting phenomenology, they 
are believed to have a minor impact upon the Universe at large, 
i.e., on the formation of structures and thermodynamics. En- 
ergetically this is not an unreasonable assumption, as VHEGR 
emission is ~' 0. 1 % of the radiative power of quasars. However, 
in this series of papers, we show that in spite of their energetic 
disadvantage, TeV blazars have a significant effect on structure 
formation and a dominant effect on thermodynamics - that is 
the VHEGR emission from blazars "punches" far above its en- 
ergetic weight. Namely, if the radiation from VHEGR sources 
is ther malized, as we have argued is the case in lBroderick et al.l 
(120121 hereafter Paper I), the heating due to these VHEGRs 
dominates photoionization heating throughout the vast major- 
ity of the Universe at z < 3, raising the temperature of the low- 
density IGM by up to two orders of magnitude. 

Given that the total power emitted by AGNs and stars in the 
UV and X-rays vastly exceeds that due to the TeV blazars, 
it seems counterintuitive that blazar heating dominates pho- 
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toheating. However, the UV and X-ray background heats 
the intergalactic medium (IGM) inefficiently after reionization, 
while VHEGR photons heats the IGM efficiently via plasma 
beam instabilities (Paper I). This difference in the heating ef- 
ficiency of photoheating vs. blazar heating is due to the dif- 
ference in the rate-limiting process in each. The rate of pho- 
toheating after reionization is not limited by the availability of 
ionizing photons, but instead by abundance of targets, and thus 
recombination. On the other hand, the heating of the IGM by 
TeV blazars does not suffer from a deficit of targets and is only 
Umited by the total cosmic power of VHEGR sources. 

We illustrate this point with the following order of magnitude 
estimate. First let us estimate the amount of photoheating that 
the IGM suffers at the present day. The recombination rate of H 
is of order the Hubble time at the mean density of the Universe 
at present. Hence, an average H atom in the IGM will recom- 
bine once over a Hubble time only to be ionized immediately 
by a UV photon. Using a spectral index of -1. 6 for the ionizin 
backg round, which is appropriate for quasars ( iFurlanetto & O, 
l2008h . the average amount of excess energy absorbed per ion- 
ization is Eexc ~ 4eV. The fraction of the rest mass energy of 
all the baryons in the universe required to produce this amount 
of heating in the IGM is 



exc ~ 



! 4 X 10" 



I / ^exc ^ 

V4eV/ 



(1) 



Hence, only a small amount of rest mass energy is injected as 
thermal energy into the IGM. An equivalent statement is that 
the diffuse IGM is optically thin to ionization radiation^ 

By comparison, the fraction of the baryon rest mass locked 
up in massive black holes in the Universe is /bh«9x 10--**and 
in stars is 0.06 (Fukugita & Peebles 2004). Assuming a radia- 
tive efficiency (relative to rest mass) for black holes and stars 
of 0.1 and 10"^ respectively, we find that fraction of the baryon 
rest mass converted to radiation in black holes and stars is 
Erad.BH ~ 10"^ and eiad.* ~ 6 X lO"*", respectivelyO This is many 
orders of magnitude larger than what is required by Equation 
([T]) and demonstrates an important point: hard ionizing radia- 
tion is inefficient at heating the IGM after H and He reioniza- 
tion. 

In Paper I, we argued that VHEGR photons are efficiently 
converted into heat in the IGM via plasma instabilities. Hence, 
the limiting factor is the total energy density of VHEGR pho- 
tons over cosmic time. To estimate this energy density, we 
note that the local observed TeV blazar luminosity density is 
2.1 xlO^^ that of the local quasar luminosity density (see Sec- 
tion [TTTT]), after correcting for various selection effects. As- 
suming that the TeV blazar luminosity density tracks the quasar 
luminosity density over cosmic time, this implies that the frac- 
tion of baryon rest mass that is converted to VHEGR photons 
(and ultimately heating of the IGM) is 

TeV Blazar Luminosity Density „ 

/xev = —p: z : 7-^ X e,ad.BH = 2.1 x 10-^ 

Quasar Lummosity Density 

(2) 

which is nearly five times that due to photoheating (i.e.. Equa- 
tion ((TJ) and demonstrates the dominance of TeV blazar heat- 

^ By "diffuse" IGM, we discount Lyman limit systems. 

^ Here we liave likely significantly overestimated the radiative efficiency 
of the stellar component for two reasons. First, most of the mass locked up 
in the stellar component is contained in low-mass stars, which are capable of 
converting only a small fraction of their rest mass into energy. Second, we have 
not accounted for the fraction of stellar radiation that is capable of ionizing the 
IGM. 



ing. Equivalently, the greater efficiency of TeV blazar heating 
compared to photoheating more than makes up for its energetic 
disadvantage. 

The physics of this heating and its cosmological conse- 
quences is the subject of this series of three papers. In Pa- 
per I, we studied the physics of VHEGR photon propaga- 
tion through the Universe. As these VHEGRs propagate 
through the Universe, they interact with the soft photons that 
comprise the extragalactic background light (EBL) and pro- 
duce ultra-relativistic pair s (see, e.g., Gould & S chreder 1 96% 
iSalamon & S tecke^ [T998t iNeronov & Semikozll2009l) . Typi- 
cal mean free paths are between 30Mpc and IGpc, depend- 
ing upon the energy of the VHEGR and redshift, i.e., the Uni- 
verse is optically thick to VHEGR. The result is a ubiqui- 
tous population of ultra-relativistic pairs, with typical Lorentz 
factors of 10^-10^. Previously, it has been assumed that 
they lose energy exclusively through inverse-Compton scat- 
tering the cosmic microwave background (CMB) and EBL, 
producing GeV gamma rays that fo r m part of the EGRB 
(see, e .g.. iNarumoto & Totanil 120061: iKneiske & MannheimI 
120081: llnoue & Totanil 120091: IVentersI 120 1 The non- 

observation of this GeV gamma-rays has been used to argue 
for cosmologically interesting IGMFs (INeronov & Vovkll20I0l: 
Tavecchio et al.ll2qTol I201II: iDermer et al.ll2011l: iTavlor et af 



201 1 ^ lDolagetal.l I2011t iTakahashi et all 1201 2tl^^wketar 
201 



We then presented a plausible alternative mechanism for ex- 
tracting the kinetic energy of the ultra-relativistic pairs: plasma 
beam instabilities. Despite the extraordinarily dilute nature of 
this ultra-relativistic pair plasma, we found a variety of plasma 
instabilities which grow on timescales short in comparison to 
the inverse-Compton cooling time, the most im portant of which 
is the "oblique" instability (iBret et al.ll2004l) . Via this insta- 
bility, these ultra-relativistic pairs lose their kinetic energy by 
depositing it as heat in the IGM. Because these beams then 
cool well before an inverse-Compton cascade (ICC) can de- 
velop, the simplest versions of the argument used to produce 
limits upon the IGMF are precluded. Hence, the existence of 
the IGMF does not follow from the non-observation of GeV 
gamma rays from existing TeV blazars as previous groups have 
argued. In addition, the lack of an ICC allows for a large and 
evolving blazar population without upsetting the Fermi limits 
on the EGRB and statistics of high-energy blazars. 

Based upon the plasma-instability mechanism, we now adopt 
as a hypothesis that the ultra-relativistic pairs primarily deposit 
their energy in the IGM via this or a related mechanism. In 
this paper (Paper II), we explore the impact of this heating on 
the thermodynamics of the IGM. We estimate the amount of 
heating provided by the observed TeV blazar population after 
correcting for the selection effects of the current pencil beam 
VHEGR observations using the all-sky monitoring of the Fermi 
satellite. We will show that the luminosity density in VHEGRs 
is of order 0.2% of the quasar luminosity density, which domi- 
nates the photoheating rate at low z- We then explore the quali- 
tative and quantitative nature of this heating, which in essence, 
serves as an alternate feedback mechanism. In particular, TeV 
blazars deposit heat evenly in a volumetric sense, i.e., indepen- 
dent of the local IGM density. Hence, this heating deposits 
more energy per baryon in low-density regions than in high- 
density regions, naturally producing an inverted temperature- 
density relation in voids. With only a minor rescaling of the 
empirically normalized of observed blazar heating, we find that 
it is possible to reproduce the inferred inverted temperature- 
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densi ty relation at z = 2 - 3 jBolton & Beckedl2009l: IViel et alj 
l2009h . something which has proven to b e a problem within the 
context of standard reio nization models dMcOuinn et al.ll2009l; 
iBolton & Becker" 2009l). 

In Pfrommer e t alJ (120 12L hereafter Paper III), we will ex- 
plore the result of this additional IGM heating upon the forma- 
tion of structure in the Universe. In particular, we will show 
that the injection of entropy into the IGM by TeV blazars con- 
tributes to developing a redshift dependent entropy floor for 
galaxy clusters and groups at z < 2 and suppresses the forma- 
tion of dwarfs. We will highlight that the redshift dependent 
nature of TeV blazar heating in our model suggests a large 
injection of entropy around z ~ 1, which boosts the entropy 
of late forming objects. This predicted enhanced entropy of 
young groups is consistent with recent observations that show 
optically bright and therefore young, groups and clusters are 
X-ray dim - that is having a lower gas density due to a raised 
entropy floor We also will show that TeV blazar heating sup- 
presses the formation of late forming d warfs both in ga lactic 
halos, i.e., the missing satellite prob lem jKravtsov 120101) . and 
in voids, i.e., the void phenomenon ( lPeeblesll2001 ) by raising 
the temperature of the IGM such that gas cannot collapse to 
form galaxies. 

This work is organized as follows: We first review the fate 
of energy carried by VHEGR photons in Section |2] discussed 
in detail in Paper I. We describe how VHEGR photons produce 
pairs in the IGM and how these ultra-relativistic pair beams 
are unstable to plasma instabilities. In particular, we highlight 
the "oblique" instability, which is especially efficient at con- 
verting the kinetic energy of the beams into thermal energy in 
the IGM. Motivated by our review of Paper I, we adopt the as- 
sumption that the kinetic energy of the ultra-relativistic pairs 
is thermalized in the IGM, either via the "oblique" instability 
or some related mechanism. In Section |3] we estimate the cur- 
rent Te V-blazar IGM heating rate by collating the known extra- 
galactic TeV blazars with a well measured spectrum, account- 
ing for incompleteness (Section [3. 1.1b . We use the similarity 
between the luminosity functions of nearby quasars and TeV 
blazars found in Paper I, to extend the heating rate to z > and 
estimate the Te V-blazar covering fraction (Sections 13.1.21 and 
I3.2l i. The implications for the thermal history of low-density 
regions (less than 10 times the mean density, i.e., 1 + (5 < 10) 
are explored in Section |4] Generally, we find that without any 
fine tuning it is possible to reproduce the inverted temperature- 
density relation at z = 2 - 3 inferred by high-redshift Lya stud- 
ies (Bolton et al. 2008; Viel et al. 200^, while simultaneously 
satisfying the te mperature constraints at z = 2 (e.g., those by 
iLidz et an201Cft and leaving the local Lya forest unaffected. 

The results of this work and Paper III assume that the energy 
of TeV blazars are efficiently thermalized in the IGM. Though 
we have identified a particularly promising instability, i.e., the 
"oblique" instability, in Paper I, the particular details by which 
the energy in ultra-relativistic pairs is thermalized are unim- 
portant. Instead, the results of this work and Paper III depend 
solely on the gross energetics of the TeV blazars. As a con- 
sequence, observations of the thermal history and constraints 
upon the structures in low-density regions represent an inde- 
pendent empirical probe of the fate of the ultra-relativistic pairs 
produced by TeV blazars. 

For all of the calculations presented below (and in this series) 
we have assumed the WMAP7 cosmology w ith ho = 0.704, 
nDM = 0.227, Ob = 0.0456, and = 0.728 (iKomatsu et al.1 
I2OIII) . 



2. REVIEW OF VHEGR PHOTON PROPAGATION THROUGH THE IGM 

The observed extragalactic VHEGR sources are all located 
at low redshift (z ^ 0.5). This is a due to the annihilation 
of VHEG Rs on the EBL, producing ultra-relativistic pairs 



. roaucing 1 

(see, e.g., 'Gould & Schreder"1967 nSalamon & Steckeril 19981: 



Neron ov & Semikoz 2009). Here we describe the fate of these 
VHEGRs, and the consequences of propagating through the in- 
tergalactic medium (IGM) for the pairs they produce. We refer 
the interested reader to Sections 2 and 3 of Paper I for a more 
complete discussion of the properties of the generated pairs and 
the importance and nature of plasma beam instabilities upon 
their propagation. 

VHEGRs are attenuated by the pair production off EBL pho- 
tons. Namely, when the energies of the VHEGRs (E) and the 
EBL photon (^ebl) exceed the rest mass energy of the pair 
in the center of mass frame, i.e., 2E Eebl(^ - cos 9) > Am^c'^, 
where 9 is the relative angle of propagation in the lab frame, 
an pair can be produc ed with Lorentz factor 7 ~ E /InieC^ 
(iGould & SchredejlI967h . An estimate for the mean free path 
of VHEGR photons is (Paper I) 



Dpp(£,z) = 35 



ITeV 



1+z 
2 



Mpc, 



(3) 



where the redshift evolution is due to that of the EBL alone, 
is dependent predominately upon the st ar formation history , 
and C = 4.5 for z < 1 and C = for z > 1 dKneiske et al.ll2004t 
iNeronov & Semikozl2009lrl . Relative to the Hubble length, the 
attenuation length of VHEGRs is very short. In fact, above 
lOOGeV the Universe is optically thick to sources at z > 1 (cf. 
iFranceschini et alj|2008l) . 

Since Dpp is much larger than any conceivable source size, 
and /iEBL ^ E, locally these pairs necessarily constitute a cold, 
highly anisotropic beam. The production of pairs by the inter- 
action of VHEGRs with the EBL is opposed by the removal of 
these pairs by various cooling processes. That is, the evolution 
of the characteristic pair-beam density at the injection Lorentz 
factor, Hb, is governed by the Boltzmann equation: 



drib c dnrib „ 
at r- or 



(4) 



where the left-hand side assumes all the pairs are moving away 
from the VHEGR source relativistically (v' = c and p'' = •yniec), 
the right-hand side corresponds to pair production, and F is the 
cooling/removal rate of these pairs. In a homogeneous steady 
state, the rate of production is balanced by the rate of removal, 
which gives: 

2Fe 

"b - j—f , (5) 

where the rate of production is given by «beam = 
2{EdN/dE)/Dpp = IFe/D^^. Generally, F is a function 
of energy and beam energy and refer the reader to Paper I for 
additional details. However, for any choice of F, the solution 
to Equation (|5]l gives n\,{E ,Fe,z)- 

Commonly, it is assumed that the pairs evolve primarily due 
to an ICC which deposits the energy from the original VHEGR 
photon near ~ lOOGeV. Note that this is fundamentally radia- 
tive; after the VHEGRs are scattered down to lOOGeV they ef- 
fectively decouple from the Universe. ICCs of these pair beams 

Despite the fact that the EBL contribution from starbursts peaks at z = 3 
and declines rapidly afterward, galaxies and Type 1 AGNs compensate for th e 
lost flux until z = 1 . See, for example. Figure 3 from lFranceschini et al.1 120081) . 



4 



CHANG, BRODERICK, & PFROMMER 



widely exploited as a possible probe of IGM magnetic fields in 
the conte xt of the missing inverse Compton fea tures at lOOGeV 
(see. e.g.jNeronov & VovB2010HT^v ecchio et al.ll201d l20Tl 

m 



Dermer et al. '20111; iTaylor et al.l 120111: iTakahashi et 



20121: 



Dolag e t al. 201 1). The associated cooling rate for this process 
for pairs with Lorentz factor 7 is 



ic ■ 



40£MCMB 

3m„c 



-7- 



1.4x 10"'°(l + z)S s"' , 



(6) 



where aj denotes the Thompson cross section. The strong red- 
shift dependence arises from the rapid increase in the CMB 
energy density with z (mcmb oc (1 +z)'^)- 

However, this may not be the dominant process. As we 
showed in Paper I, plasma beam instabilities are a potential 
mechanism by which the kinetic energy of the pairs is extracted 
much more rapidly. The two plasma beam instabilities most 
frequently discussed, the two stream and Weibel instabilities, 
are particular limits of a more general instability, differentiated 
by the direction of the perturbed wave vector with respect to 
the beam or ientation (perpen dicular for Weibel, parallel for the 
two stream: lBret et al.ll2005h . For dilute beams, by far the most 
powerful growth occurs at oblique angles, and thus na med the 
"oblique" instability. Paper I and iBret et all (120101^ gives a 
more extensive discussion of these three plasma instabilities. 

Relative to the laboratory frame, these pair beams are "cold", 
i.e., their transverse momentum is much smaller than their par- 
allel momentum, but even the small transverse temperatures 
that are acquired from pair production for these pair beams are 
important, placing the oblique instability in the kinetic regime. 
Here, the oblique instabiUty cooling rate has been numerically 
measured to be 



M,k : 



:0.4- 



«b 



kTh niGM 



Ldp ~ 0.47-1 



/ 47re-«g 



(7) 



where nicivi = 2.2 x 10"^(1 +(5)(1 +z)^cm"^ is the IGM free- 
electron number density (assuming full ionization), e is the el- 
ementary charge, and we have set the beam temperature in the 
beam frame to nieC^/k, characteristic of that induced by pair 
production, and thus kTt, = mpC^/7 jBret et aI1l2010ah . Both 
the cold and hot growth rates have been verified explicitly us- 
ing particle-in-cell (PIC) simulations, though at som ewhat less 
dilute beams than the pair beams from TeV blazars jBret et al.l 
l20T0bl) . 

For the kinetic oblique instability, the cooling rate is 



rM,k^3.6x 10"" (1+5) 



-1/4 



(6C-3)/4 



111 

2 

^1045 erg s-i) [jeV ) 



f ELe 



1/2 



f E \ 



3/2 



s-', (8) 



where Le is the isotropic-equivalent luminosity per unit energy 
of the VHEGR source. This is a stronger function of photon en- 
ergy than inverse-Compton cooling, implying that it will even- 
tually dominate at sufficiently high energies, assuming a flat 
VHEGR spectrum. In addition it is a very weak function of 5, 
being only marginally faster in lower-density regions, and thus 
the cooling of the pairs is largely independent of the properties 
of the background IGM. 

The effective cooling rates induced by the "oblique" insta- 
bility can be obtained by numerically solving Equation (|5]l for 
nh, with r = Tic + rM.k as described in Paper I. For the purpose 



of this paper, for TeV blazars the effective cooling rate is gen- 
erally dominated by the "oblique" mode and thus the fraction 
of VHEGR emission that is effectively thermalized is nearly 
unity. Therefore, we are justified in assuming that a ll of the 
VHEGR emission is thermalized. In Section 13.1.21 we will 
briefly revisit this assumption at high z, finding again that it 
is well justified at all z0 

However, this may not be the case for all potential VHEGR 
sources. The physics of the "oblique" instability, like all plasma 
beam instabilities, depends strongly upon the beam density. 
For the pair beams that result from VHEGR-EBL photon in- 
teractions, the beam density is a function of source luminosity. 
Hence for sufficiently low luminosity systems, the "oblique" 
mode grows so slowly that inverse Compton off of the CMB 
dominates it. This defines a critical isotropic-equivalent lumi- 
nosity for plasma beam instabilities to be relevant, typically 
near 10^^ erg s"' , though depending upon redshift (see Figure 3 
of Paper I). For the TeV blazars considered in this paper, this is 
generally not a concern. However, this will be important in our 
discussion of alternative VHEGR sources in Section [374l 

3. TEV BLAZAR HEATING OF THE IGM 

As we have argued in Paper I and reviewed in Section |2l 
plasma instabilities on the pair beam dissipate the bulk of 
the VHEGR emission. In particular, the "oblique" instabil- 
ity appears to be a promising mechanism by which the ultra- 
relativistic pair beams that are produced from the interaction of 
VHEGR and EBL photons blazars are efficiently thermalized 
and converted to local IGM heating. Henceforth, we will make 
the assumption for the rest of this work and again in Paper III 
that the energy of TeV blazars is efficiently thermalized in the 
IGM. With this stated assumption, we now discuss the sources 
that dominate the IGM heating rate, estimate the magnitude of 
the heating rate and its evolution with redshift, and briefly as- 
sess the homogeneity with which this new process occurs. 

3. 1. TeV-Blazar Heating Rates 

The local heating rate associated with the dissipation of the 
high-energy emission from a single blazar is determined by 
two factors: the rate at which the VHEGRs are converted into 
pairs and the rate at which the energy of the pairs is subse- 
quently converted locally into heat. The former is determined 
by the pair-production cross-section. Due to the efficiency of 
the plasma instabilities in dissipating the beam energy, the lat- 
ter is set by the ratio of the relevant cooling rates. That is, the 
single-blazar heating rate is 



q= dE 



OjE) 
Dr,JE,z) 



f{FE,E,z)FE, 



(9) 



where 9{E) is a dimensionless function due to the pair-creation 
threshold and which depends upon the shape of the EBL spec- 
trum (we set 9{E) to vanish for E < lOOGeV and be unity oth- 
erwise), and 



f(FE,E,z)=l-fi 



M.k 



IC ■ 



Tic + Tjyi.k 



(10) 



which is a function of Fe, E, and z via the dependence of the 
cooling rates upon «b, 7, and z. Generically, / represents the 

' We, however, have made the implicit assumption that the nonhnear state 
of the "oblique" instability removes kinetic energy from the beam at the growth 
rate of the linear instability. 
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fraction of pair energy that is thermalized, but we have chosen a 
specific form of / corresponding to the "oblique" instability to 
make the discussion below more concrete. While it is also very 
weakly dependent upon 6, this may be neglected in practice. 

Within the linear regime the plasma instabilities responsible 
for the dissipation of pair beam energy are independent of those 
arising from beams in substantially different directions. Since 
the local VHEGR flux is dominated by a number of sources that 
is small in comparison to that needed to isotropize the beam's 
phase space, we may treat the resulting evolutions of their asso- 
ciated pair beams independently and sum their resulting heat- 
ing rates to determine the total heating rate, Q. 

3.1.1. Estimating the Local Heating Rate 

The local heating rate can be estimated in at least two ways, 
both of which give similar results. First, since the high- 
energy gamma rays deposit their energy locally, we can iden- 
tify the local heating rate with the high-energy gamma-ray 
luminosity density of TeV blazars. In Paper I we showed 
that this is roughly 2.1 x 10"^ that of quasars, and thus 
approximately (0.5-1.4) x lO^'^ergMpc"^ s"' = (3.4-9.7) x 
10~^eVcm"^Gyr"', depending upon the minimum luminosity 
at which the heating mechanism operates. 

Alternatively, given a sufficiently complete sample, we can 
estimate the local heating rate using that implied by the fluxes 
of the observed TeV blazars. In the present epoch, /{Fe ,E,z) — 
1 at the relevant E, and thus the heating rate takes the particu- 
larly simple form; 



(11) 



To evaluate this sum, we have collated the presently 46 extra- 
galactic VHEGR sources knownQ Of these 46, only 28 have 
published measurements of their VHEGR flux from a combi- 
nation of VERITAS, H.E.S.S., and MAGIC observations. For 
these 28 sources, we have extracted the parametrized spectra 
assuming the form. 



dN 
'dE 



■■fo 



Eo 



(12) 



where fo is the normalization in units of cm' -^s-'TeV"'. The 
gamma-ray energy flux is trivially related to dN /dE by Fe = 
EdN /dE cx £■'"", from which we obtain a VHEGR flux. 



F = Eofo 



lOTeV 



100 GeV 



dE ( — 

Eo 



l-a 



(13) 



and for sources with a measured redshift a corresponding 
isotropic-equivalent luminosity, L = AnDj^F, where Dl is the 
luminosity distance. The resulting fo, Eo, a, F, and L are col- 
lected in Table [T] In addition we list the redshift, inferred dis- 
tance, and absorption-corrected intrinsic spectral index at Eo, 

* We should note that the if the number of blazars in the sky is so numerous 
as to make the phase space distribution of beam particles isotropic, certain 
classes of plasma instabilities - in particular the Weibel instability - will be 
suppressed. The volume filling factor of a pair beam in phase space is the ratio 
of perpendicular beam temperature to the beam energy. kT\,/(me(r) ~ 10"* 
(see Paper 1); hence we would need Af ~ lO'" blaz ars to isotropize any given 
point in space. However, as will be shown in Section [331 the number of blazars 
that contribute substantially to the local heating rate is much smaller; hence we 
believe the local pair distribution function is sufficiently anisotropic. 

^ See |http://www.mppmu.mpg.de/~rwagner/sources/ | for an up-to-date list. 
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Figure 1. Local heating rate due to only the observed TeV blazars as a function 
of upper-cutoff (fimax) and lower-cutoff energies (red, black, and blue lines). 
In particular, note that the overall value of gobs is relatively insensitive to rea- 
sonable variations in the spectral cutoffs. 



obtained via 



dlnE 



-t(Eo,z) , 



(14) 



Eo 



where t(E,z) is the optical depth accrued by a VHEGR emitted 
at redshift z and with observed energy E^ For high-redshift 
sources a can be less than 2, implying that an intrinsic spectral 
upper-cutoff must exist. Here we conservatively take this to 
be at £ ~ lOTeV, which is well justified given the distances 
to the two sources that dominate the observed local TeV flux, 
Mkn 421 and lES 1959H-650, though as we shall see below, the 
heating rate is relatively insensitive to the particular values of 
the lower and upper spectral cutoffs. 

Using the observed VHEGR source in Table [T] the heating 
rate associated with the 23 IBL/HBL blazars (labeled I and H 
in the table, see below) is 



obs - 



obs AGN 



Eofo 



10 TeV 



dE 



Dpp{Eo,0) JiooGeV ^0 



-V 

EoJ 



(15) 



Note that for a ~ 3, characteristic of the two dominant TeV 
sources, this depends logarithmically upon the lower and upper 
spectral cutoffs. The resulting gobs is shown in Figure [1] as 
a function of the upper spectral cutoff, E,^^^ for a handful of 
different lower cutoffs, ranging from 50 GeV to 200 GeV. For 
our fiducial values, Qobs = 1.2 x 10~^eVcm"^Gyr"', with this 
only weakly depending upon the various cutoffs, changing by 
at most by a factor of 2 over the reasonable ranges. Hence, for 
the remainder of this work, we choose a lower spectral cutoff 
of 100 GeV. 

To complete this estimate, we now correct for the selection 
effects of VHEGR observations. To correct for the pointed 
nature of VHEGR observations, we rely on the all- sky GeV 
gamm a-ray observations from the Fermi satellite (lAbdo et al.l 
l2010al) of TeV blazar. Those belong to 2 subclasses of blazars, 
namely high-energy peaked BL Lacs (HBL) and the somewhat 
less efficient accelerators, intermediate-energy peaked BL Lacs 

* This differs subtly from the definition of te(E,z) in Paper 1, where there 
we set E to the emitted energy of the gam ma ra y. A full definition at arbitrary 
observer redshift can be found in Equation {A2\ . 
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Table 1 

List of TeV Sources with Measured Spectral Properties in Decreasing lOOGeV-lOTeV Flux Order 



Name 














log 10 


a s 




Class 


Mkn421 


0.030 


129 


68 


1 


3.32 


1.7 X 10^ 


45.6 


3.15 


44 


H 


lES 1959+650 


0.047 


201 


78 


1 


3.18 


1.6 X 10' 


45.9 


2.90 


47 


H 


lES 2344+514 


0.044 


190 


120 


0.5 


2.95 


2.3 X 10^ 


45.0 


2.82 


9.2 


H 


Mkn 501 


0.034 


150 


8.7 


1 


2.58 


85 


44.4 


2.39 


6.0 


H 


3C 279 


0.536 


2000 


520 


0.2 


4.11 


68 


46.9 


2.53 


1.0 


Q 


PKS 2155-304 


0.116 


490 


1.81 


1 


3.53 


64 


45.4 


2.75 


1.4 


H 


PG 1553+113 


> 0.09 


> 380 


46.8 


0.3 


4.46 


41 


>44.9 


<4.29 


< 3.88 


H 


W Comae 


0.102 


430 


20 


0.4 


3.68 


31 


44.9 


3.41 


0.6 


1 


3C 66A 


0.444 


1700 


40 


0.3 


4.1 


28 


46.3 


2.43 


0.4 


1 


lES 1011+496 


0.212 


870 


200 


0.2 


4 


26 


45.5 


3.66 


0.4 


H 


lES 1218+304 


0.182 


750 


11.5 


0.5 


3.07 


24 


45.4 


2.37 


0.8 


H 


Mkn 180 


0.045 


190 


45 


0.3 


3.25 


20 


44.0 


3.17 


0.6 


H 


IH 1426+428 


0.129 


540 


2 


1 


2.6 


20 


45.0 


1.71 


1.4 


H 


RGB J07 10+591 


0.125 


520 


1.36 


1 


2.69 


15 


44.8 


1.83 


0.9 


H 


lES 0806+524 


0.138 


580 


6.8 


0.4 


3.6 


10 


44.7 


3.21 


0.2 


H 


RGB J0152+017 


0.080 


340 


0.57 


1 


2.95 


8.5 


44.1 


2.45 


0.3 


H 


lES 1101-232 


0.186 


770 


0.56 


1 


2.94 


8.2 


44.9 


1.50 


0.3 


H 


lES 0347-121 


0.185 


770 


0.45 


1 


3.1 


8.2 


44.9 


1.67 


0.3 


H 


IC 310 


0.019 


83 


1.1 


1 


2.0 


8.1 


42.8 


1.90 


0.1 


H 


PKS 2005-489 


0.071 


300 


0.1 


1 


4.0 


8.0 


44.0 


3.56 


0.1 


H 


MAGIC J0223+430 






17.4 


0.3 


3.1 


7.6 




< 3.1 


0.2 


R 


lES 0229+200 


0.140 


590 


0.7 


1 


2.5 


6.4 


44.5 


1.51 


0.5 


H 


PKS 1424+240 


< 0.66 


< 2400 


51 


0.2 


3.8 


6.3 


<46.1 


> 1.42 


0.1 


1 


M87 


0.0044 


19 


0.74 


1 


2.31 


5.9 


41.4 


2.29 


0.6 


R 


BL Lacertae 


0.069 


290 


0.3 


1 


3.09 


5.4 


43.8 


2.67 


0.2 


L 


H 2356-309 


0.165 


690 


0.3 


1 


3.09 


5.4 


44.6 


1.86 


0.2 


H 


PKS 0548-322 


0.069 


290 


0.3 


1 


2.86 


4.0 


43.7 


2.44 


0.2 


H 


Centaurus A 


0.0028 


12 


0.245 


1 


2.73 


2.8 


40.7 


2.72 


0.2 


R 



Reference 



Chandra et al. (2010) 

Aharonian et al. (2003) 

Albert et al. (2007c) 

Huang et al. (2009) 

MAGIC Collaboration et al. (2008) 

HESS Collaboration et ak (2010) 

Aharonian et al. (2008b) 

Acciari et al. (2009b) 

Acciari et al. (2 009c) 

Albert et al. (2007b) 

Acciari et al. (2010a) 

Albert et al. (2006) 

Aharonian et al. (2002) 

Acciari et al. (201()c) 

Acciari et al. (2009a) 

Aharonian et al. (2008a) 

Aharonian et al. (2007a) 

Aharonian et al. (200713) 

Aleksic et al. (2010) 

Aharonian et al. (2005) 

Aliu et al. (2009) 

Aharonian et al. (2007^ 

Acciari et al. (2010b) 

Acciari et al. (200|) 

Albert et al. (2007a) 

Aharonian et al. (2606K) 

Aharonian et al. (2010) 

Raue et al. (2010) 



Comoving distance in units of Mpc 

Normalization of the observed photon spectrum thai wc assume to be of the form dN jdE - fQ{E j Eq)~'-^ , in units of 10~ 
Energy at which we normalize the spectrum, in units of TeV 
^ Observed spectral index at Eq 

^ Integrated energy flux between 100 GcV and 10 TeV. in units of 10"^" ergcm"-s"^ 
^ Inferred isotropic integrated luminosity between 100 GeV and lOTeV, in units of ergs~^ 
^ Inferred intrinsic spectral index at Eq 

^ Local plasma-instability heating rate in units of 10~"^cVcm~^Gyr~^ 

^ H, I, L, Q, and R correspond to high-energy, intennedi ate -energy, low-energy peiiked BL Lacs, flat spectrum radio quasjirs, 

(IBL) which are, in some cases, also able to reach energies 
beyond 100 GeV. Outside of the Galactic plane, Fermi ob- 
serves 118 high-synchrotron peaked (HSP) blazars and a total 
of 46 high-synchrotron peaked (ISP) blazars0 Roughly half 
of the latter are likely to emit VHEGRs as indicated by their 
flat spectral index between 0.1 and 100 GeV, F < 2 (see th e 
spectral index distribution of Figure 14 in lAbdo et al.l (l2009al) ). 
Of these potential 141 TeV blazars, only 22 have also been 
coincidentally identified as TeV blazars (out of a total of 28 
coincident TeV sources), while there are a total of 33 known 
TeV blazars (29 HBL, 4 IBL). If these 141 sources are all TeV 
emitters, but have not been detected due to incomplete sky cov- 
erage of current TeV instruments, then the selection factor is 
??sei = 141 /33 = 4.3. In addition, the duty cycle of coincident 
GeV and TeV emission is r/duty = 33 /22 = 1.5. Here we assume 
that the luminosity distribution of observed VHEGR sources 
reflects the true distribution after correcting for the effects of 
flux limitations in the observations. That this may be done is 
justified empirically in Section 5.1.2 of Paper I, and demon- 
strated explicitly in Figure 5 of Paper I. Thus, we may use 
constant correction factors (independent of luminosity) to es- 
timate the true distribution. Finally, by excluding the Galactic 
plane for galactic latitudes b < 10°, this is an underestimate by 
roughly Tysky = 1.17. Taken together, the true heating rate in our 

' The source classes of HSP/ISP are very similar to the commonly used 
HBL/IBL classes. Hence we identify both for the remainder of this work. 



and radio galaxios of Fjaanoff-Riky Typo I (FR I), rospoolivoly. 

"standard" model is then: 

e|,_„~7xio-«(^) 

where rj^y^ is a remaining coefficient of order unity correcting 
for systematic uncertainties. These include corrections to the 
spectral model we currently use and the corrections to the com- 
pleteness of the Fermi sample of HBLs and IBLs in accounting 
for the complete population of VHEGR sources. The fact that 
there are already 4 radio galaxies of Faranoff-Riley Type I (FR 
I), 2 flat spectrum radio quasars, and 2 starburst galaxies emit- 
ting VHEGRs implies that we are probably too conservative 
in accounting for all the VHEGR sources. To explore the ef- 
fects of the variation in the heating rate, we also adopt an "op- 
timistic" heating model that is normalized t o fit the observe d 
IGM inverted temperature-density relation (IViel et al.ll2009h . 
This "optimistic" model has 77sys = 1.6 yielding a heating rate 
of G|_=o = 1-4 X 10-^eVcm-^^Gyr-i. 

In estimating Q from Equation (fTTl i we have made a num- 
ber of implicit assumptions. First, we have assumed that the 
observed distribution of TeV blazars locally is representative 
of the average distribution at any given point in the Universe. 
Since the blazars that have been observed locally go out to 
z ~ 0.5, we believe this is a relatively safe assumption. Sec- 
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ond, we have assumed that the current sample of blazars is suf- 
ficiently flux-complete to dominate the heating rate at Earth. 
That is, we have assumed that the effective flux limit of the 
current generation of Cerenkov telescopes is low enough that it 
captures the bulk of the sources responsible for the local heat- 
ing of the IGM. T hat th is is the case is less clear, and will be the 
subject of Section[32]and the Appendix in some detail. Unfor- 
tunately, since the constant of proportionality relating the local 
TeV blazar and quasar luminosity densities is obtained using 
the observed set of blazars, our two estimates are strongly cor- 
related. 

Nevertheless, this provides a convenient lower limit upon Q, 
and we use this heating rate, which we denote the "standard" 
heating rate, as the mean heating rate of any fluid element in 
the present-day Universe. Over a Hubble time, the total heat 
deposited per unit volume is then roughly 



eL=, 



Mtot 



Z=0 



;=0 



10"''eVcm"^ 



(17) 



which is sufficient to raise the temperature of 1 + (5 = 0. 1 regions 
to approximately 1.7 x lO^K, and thus dramatically alter the 
thermal history of voids. 

3.1.2. Estimating the z > Heating Rate 

Extrapolating the above estimate to z > requires under- 
standing how the TeV blazar properties and population have 
evolved, as well as the evolution in f(FE,E,z). Generally, the 
average heating rate is given by 

f ~ n 

Q= dVdlogiQLdadn(l)B{z;L,a,n)—q, (18) 
J 27r 

where (/)B(z;L,a,ri,z) is the physical number density of blazars 
at a given redshift per unit logarithmic isotropic-equivalent lu- 
minosity, spectral index, and blazar jet opening angle (we've 
assumed all jets are symmetric). We make the simplifying 
assumption that (/>£ is separable into components describing 
the evolving luminosity density distribution, (j)B(z,L), and a 
static, unit-normalized spectral distribution, (pB{a,^l), where 
(^B(z;L,a,f2) = (j)B{z,L)(pB{a,n). This produces. 



Q= / c/logi()LL(/)B(z,L) 



X / D^dDdadn / dE 



where the bulk of the redshift and luminosity dependence is 
now contained in the outer-most integral over logjgL, and we 
have used rg ~ D/Dpp (since the heating is dominated by 
nearby objects, for simplicity we do not distinguish between 
different cosmological distance definitions, though see the Ap- 
pendix for a more careful treatment). If there were no redshift 
or flux dependence in the remaining terms, it would be possible 
to simply normalize the heating rate by Q | _^ and the estimated 
TeV-blazar luminosity density, which we define to be 



Ab(z) ; 



dlogiQLL(j>B(z,L), 



(20) 



log,oi-m 



where Ln^n — 3 x lO'^^ergs ' is chosen so that the plasma insta- 
bilities operate efficiently at all redshifts of interest. However, 




l+z 



Figure 2. Heating correction due to the intrinsic evolution of tiie plasma in- 
stability and inverse-Compton cooling rates as a function of redshift. Shown 
are the corrections for L = lO^'^ergs"' (black solid), lO^'^ergs"' (blue long- 
dashed), and 10'*'' ergs"' (red short-dashed). For comparison, the TeV blazars 
that dominate the heating locally have luminosities of roughly 5 x 10*^ ergi"'. 

this is not entirely the case since f{FE,E,z) retains some de- 
pendence upon z, and thus we set 



^-A,(0)^l-«®^ 



Az), 



(21) 



where Qcorr provides a correction due to the changing strength 
of the pair beam dissipation mechanism in comparison to 
inverse-Compton cooling. 

The VHEGR flux at Earth is dominated by a handful of 
very bright sources for which the absorption corrected a ~ 3 
at ITeV and isotropic-equivalent luminosity 5 x Iff^^ergs"'. 
This simplifies the estimation of Qcorr significantly, giving 



Scon-(Z) : 



dEdD- 



OE- 



-D/D„ 



pp 



,E,z 



47rD2 

dEOE-^ , (22) 



in which the form and magnitude of Le is fixed. The 
resulting correction factors are shown in Figure |2] for 
isotropic-equivalent luminosities above lOOGeV ranging from 
lO'^'^ergs"' to 10"^^ ergs"'. Generally, Qcon makes a small 
(< 15%) correction following the peak of the quasar luminosity 
function {z = 2), though can grow to as much as 80% by z = 6 
for dim objects. However, since we will be primarily interested 
in bright blazars at z < 4, Qcorr modifies the heating rate at rel- 
evant redshifts by < 50% (see also Section [33] l. Thus, given 
the comparatively larger uncertainties associated with the TeV 
blazar luminosity function, plasma heating mechanism, and ad- 
ditional potential VHEGR sources, we neglect Qcon in what 
follows. 

Therefore, obtaining the heating rate as a function of z is 
reduced to estimating Ab(z)/Ab{0), the normalized evolution 
of the blazar luminosity density. However, no VHEGR emit- 
ting blazars are known beyond z — 0.7, presumably due to the 
pair-creation absorption associated with propagation through 
the EBL. As a consequence, nothing is known about Ab(z) at 
high redshifts directly. Nevertheless, in Paper I we showed that 
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Figure 3. Blazar luminosity density (in comoving units, with Tjsys = 1) as 
a function of redsliift. The shaded region represents the l-cr uncertainty that 
results from a combination of the uncertainty in the number of bright blazars 
that co ntribute to the local heating and in the uncertainties in the Hopkins et al] 
<2007l) quasar luminosity density to which we normalize. Our optimistic and 
standard models, defined in Section [3.1.1 1 are shown by the magenta short- 
dash-dot and blue long-dash-dot lines. Insert: Comparison between the ob- 
served blazar luminosity function (data points with one sigma eiTor bars and 
upper limits) with the quasar luminosity function that has been shifted by a 
factor of 2.1 X 10"' in luminosity density and 0.55 in luminosity. This insert 
is a simplified reproduction of Figure 5 of Paper I. 



a close relationship between and the quasar luminosity func- 
tion of iHopkins et al.l ( 120071) . 0q, exists at low z: 



/)b(0.1,L)~3.8x 10-^(?iQ(0.1,1.8L). 



(23) 



The insert in Figure |3] shows a comparison between the ob- 
served TeV-blazar luminosity function, calculated from the dis- 
tribution of the sources in Table [T] and applying an empiri- 
cally determined flux limit of 4.19 x 10~'^ergcm~^s~', with 
(/)b(0.1,L) (for a full description of how the luminosity func- 
tion was formed, and how it compares to (^q, see Section 5 of 
Paper I). Furthermore, we showed that once the ICCs are sup- 
pressed by the plasma beam instabilities (or some analogous 
mechanism), extending this relationship to high-z was in excel- 
lent agreement with the best current constraints upon the high- 
z TeV blazar population: the Fermi TeV blazar statistics and 
the Fermi measurement of the extragalactic gamma-ray back- 
ground (between lOOMeV and lOOGeV). Thus, we estimate 
Ab(z) — 2.1 X 10~-'Ag(z), shown in Figure |3] where Aq(z) is 
the luminosity density of quasars. 

Figure|3]shows that employing the quasar luminosity density 
to estimate the heating rates has profound consequences. In 
particular, the inferred luminosity density of TeV blazars 
(solid line) rises rapidly with increasing z, with the comoving 
density increasing by r oughly a factor of ^ 10 by z = 2 
(see also Figure 8 of IHopkins et al.l 12007b . In physical 
units this corresponds to a increase by a factor of nearly 
300. Thus, we expect an increase in the local heating rate 
by a similar factor over the presently estimated value near z = 2. 
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Figure 4. Comoving mean separation of blazars with isotropic-equivalent 
luminosities (L, in the lOOGeV-lOTeV band) above L,„ = lO^'^ergs"' (dot- 
dash), 10'*'' ergs"' (solid), 10''^ ergs"' (short-dash), 10''* ergs"' (long-dash), 
and the local luminosity-weighted median luminosity of TeV blazars (red dot), 
as functions of redshift. For reference the spectrally averaged (with a = 3) 
mean free path, Dpp, and the local mean free path for a 1 TeV gamma ray are 
shown by the grey solid and dashed lines, respectively. For sources with a 
spectral break between lOOGeV and ITeV, we anticipate the effective mean 
free path to lie between these. 



3.2. Homogeneity of TeV-Blazar Heating 

We now investigate the assumption of even heating as a pre- 
lude to our simple one-zone model of the IGM. In lieu of large- 
scale simulations, the homogeneity of the heating due to TeV 
blazars is difficult to assess for a variety of reasons. First, the 
duty cycle of the TeV blazars is unknown. Second, the density 
at high redshifts is poorly constrained. Third, the importance of 
clustering bias is unclear. Fourth, it is difficult even to define 
which blazars are relevant, e.g., which luminosity range con- 
tributes the bulk of the local heating. Nevertheless, we make 
an attempt to roughly characterize the possible range in the 
stochasticity of the local heating rates via a number of different 
estimates. 

3.2.1. Mean Separation of TeV Blazars 

We begin by estimating the instantaneous comoving mean 
separation of visible blazars above various luminosity thresh- 
olds. 



Db(L„„z) = (1+z) 



log]„ Lm 



-1/3 



(24) 



where the precise value of Lm is unimportant as long as it is 
above the peak of the luminosity function; here we choose 
Lm = 2 X lO'^^ergs"', consistent with the theoretically expected 
upper limit of TeV blazar luminosities. Since we have deter- 
mined (j)B from the observed blazar population, and are making 
use of the isotropic-equivalent luminosities, this is independent 
of the blazar jet opening angle; smaller opening angles will re- 
sult in a correspondingly larger number of objects such that the 
total number seen is unchanged, and therefore Db is fixed. A 
comparison between Db and the mean free path of VHEGRs 
is shown in Figure |4] Because Dpp varies dramatically from 
lOOGeV to 10 TeV, for this purpose, we define a spectrally av- 
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eraged mean free path, Dpp, determined implicitly by 

nlOTeV /.lOTeV 

/ ^££i-"e-flppfe)/Opp(£,z) = e-i / dEE^-°', (25) 

JlOOGeV JlOOGeV 

in which we choose a = 3 based upon the local TeV blazar sam- 
ple. This is roughly the e-folding distance of the entire spectral 
band, shown in Figure|4]by the grey solid line, and in practice is 
quite close to Dpp(225 GeV, z). If the intrinsic TeV blazar spec- 
tra peak above lOOGeV, our estimate of the typical VHEGR 
mean free path could be significantly too large; thus we also 
compare Dg to Dpp at 1 TeV (the grey dashed line in Figure 
IDi, providing an extreme lower-bound upon the spectrally av- 
eraged mean free path in practice. The rapid increase of the 
density of EBL photons with z, associated with the larger star 
formation rate in the recent past, results in a substantially re- 
duced Dpp by z = 1. Prior to z = 1, the physical number density 
of EBL photons remains nearly constant, and thus Dpp oc (1 +z) 
in comoving units. 

In the present epoch, Dpp is quite large, and thus de- 
spite their sparsity the mean separation of even bright blazars 
(10"^^ ergs"') is is less than Dpp. This remains true for ob- 
jects with luminosities < 10"*^ ergs"', which includes the lo- 
cal luminosity-weighted median TeV blazar luminosity, Lo.5(z), 
defined such that 

/•logic i-M 

/ L^B(L,z)aflogioL = 0.5AB(z). (26) 

Jlog,|,Lo.5fe) 

Hence, locally, we expect blazar heating to be quite uniform. 

At high redshift matters change somewhat. Until z = 1, 
Dppiz) and Dpp(lTeV,z) both decrease more rapidly than Db- 
For z > 1, in comoving units Z)pp(z) and Z)pp(l TeV,z) increase 
slowly, though at a marginally larger rate than Db- As a re- 
sult, near z ~ 1 the mean separation of TeV blazars is largest 
in comparison to the VHEGR mean free path. Thus, the mean 
separation between objects more luminous than 10'*^ ergs"' is 
larger than Dpp near z ~ 1 . Nonetheless, when the lower lu- 
minosity limit is dropped to 10'*^ ergs"' or below, generally 
Db(z) < Dppiz) for all z < 6. Similarly, the mean separation be- 
tween blazars more luminous than Lo.5(z) is also smaller than 
our estimate of Dpp for the relevant redshift range. Hence we 
may expect that nearly all patches of the Universe will be illu- 
minated by at least one luminous TeV blazar following z ~ 6. 

3.2.2. Estimates of the Number of TeV Blazars that Contribute 
Significantly to the Heating Rate 

Closely related to the mean separation is the number of 
blazars within the f ~ D/Dpp = 1 surface (defined explicitly in 
the Appendix) above some luminosity limit. Figure |5] shows 
these for objects with luminosities above 10^^ ergs"' (blue 
short-dash), and Lq 5 (red dot). While these may be roughly 
inferred from the associated mean separations in FigureH] they 
differ slightly from {4-n/3){Dpp/DB)^ due to the rapidly evolv- 
ing Dpp and blazar population combined with the finite look- 
back time in the integral. At z = both are well above unity, 
with the luminosity-weighted median value substantially ex- 
ceeding our estimate of roughly 170 visible extragalactic TeV 
blazarsPl However, since generally Dpp{z) — 4.4Dpp(l TeV,z) 

Recall that in Table [T] we have listed 23 TeV Blazars, but a correction 
factor of 7.5 needs to be applied to account for the incompleteness of present 
TeV surveys. 
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Figure 5. Number of TeV blazars, estimated in various ways, that contribute 
to the heating of a given patch as a function of redshift. Definitions for A/g 
include: the number of blazars with intrinsic isotropic-equivalent luminosities 
above 10'*^^ergs"' within the f = 1 surface (blue short-dash), the number of 
blazars with luminosities above that luminosity-weighted median value at each 
redshift within the f = 1 surface (red dot), and the number of blazais with 
individual heating rates that exceed that above which 50% and 75% of the total 
heating is produced (black solid and long-dash, respectively). For comparison, 
our estimate of the completeness-corrected number of TeV blazars that are 
presently observable in the TeV is shown by the filled back point, with error 
bars denoting the Poisson uncertainty only. 

not all of the blazars that contribute to the local heating are 
expected to appear as strong TeV sources. 

The low-z behavior of the number of blazars is dictated by 
the rapid evolution in Dpp, associated with the recent rapid vari- 
ation in the star formation rate (and thus the number density 
of EBL photons). Prior to z = 1, the EBL density, and hence 
Dpp, is roughly constant in physical units, and the evolution 
number of visible blazars becomes indicative of the underlying 
evolution of blazar population. At all z < 6 these estimates of 
the number of blazars responsible for the bulk of the heating, 
A/fi, exceed unity, implying only a small fractional spatial vari- 
ation in the blazar heating rate. However, while the luminosity- 
weighted median estimate of Ab does give some idea of which 
population is responsible for most of the TeV blazar luminosity 
density, neither encapsulates the population responsible for the 
majority of the local heating. Thus, if the local heating were 
dominated by a handful of luminous sources, it is possible for 
the stochasticity to be much larger 

Therefore, we show a third estimate of Ab, corresponding 
to the number of sources with ^s larger than the heating-rate- 
weighted median value (see the Appendix for a precise defi- 
nition). That is, at a given redshift, the Mb sources with the 
highest ^s generate half of the total heating rate. Thus, if a sin- 
gle source were to dominate the local heating rate, JVb would 
be considerably less than unity, indicating a large degree of 
variability. However, below z ~ 3.5 this is not the case; the 
fractional heating rate-defined Mb is considerably larger than 
unity. At high redshifts the heating rate becomes increasingly 
dominated by more distant, luminous, and rarer objects, de- 
parting from the previous estimate of Mb- 

Nevertheless, while forz ^3.5 few sources contribute nearly 
half of the heating rate, the total heating rate must be compa- 
rable to the total VHEGR luminosity density of the blazars. 
Therefore, the total number of contributing objects must be 
similar to the number of visible blazars above the luminosity- 
weighted median L (i.e., the red dotted line in Figure|5]). As a 
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consequence, the number of relevant sources must rapidly rise 
with decreasing fractional heating rate. This is explicitly indi- 
cated by the dashed line in Figure |5] which shows the number 
of sources responsible for 75% of the local heating rate. 

In any case, it is clear that for the redshifts at which blazar 
heating is likely to be important, z < 4, the heating rate will 
be relatively uniform. Between z ^ 4 it may experience or- 
der 50% fluctuations, and by z ~ 6 will exhibit order unity 
deviation. However, we note that all of our estimates of the 
heating inhomogeneity are predicated upon the assumption that 
the (f)B accurately represents the distribution of blazars, and has 
the associated considerable uncertainties. For a more thorough 
discussion of different estimates of the number of contributing 
blazars, and a discussion of which blazars dominate the heating 
rate at a given redshift, see the Appendix. 

In addition to the statistics of blazars, the homogeneity also 
depends upon the intrinsic variability of TeV blazars. During 
the time that the TeV blazars have been observed their VHEGR 
emission has remained remarkably stable. However, this pro- 
vides only a weak lower limit 4yr) upon the variability 
timescale in these objects ( iDermer et al.ll20in) . In addition to 
depending upon the properties of the source itself, the variabil- 
ity of TeV blazars depends upon the primary emission mech- 
anism responsible for the VHEGR component of TeV blazars. 
There are two classes of inverse-Compton models to explain 
the VHEGR emission: 

1. Synchrotron self-Compton (SSC) model: In this model, 
the synchrotron radiation field is Compton up-scattered 
to TeV energies by a relativistic electron po pulation 
(iJones et alJll974tlGhisellini & Maraschilll989h . Recent 
work has led to the conclusion that a simple homoge- 
neous, one-zone, SSC model cannot explain t he SEP 
of the majority of blazars (see Figure 36 of Ab do et all 
1201 Obi) . However, models with multiple SSC compo- 
nents are consistent with the blazar SEDs. Typically 
these invoke a steady, primary SSC component which 
peaks at the IR/optical (S) and 7-ray band (IC), and a 
second more energetic and usually more variable com- 
ponent, which peaks in the UV or X-ray band (S) and at 
GeV/TeV energies (IC). The variability of this energetic 
component increases the probability that a given patch of 
the IGM will see a TeV blazar during its history, result- 
ing in larger homogeneity in the blazar heating. 

2. External radiation Compton (ERC) scenario: this model 
proposes that the relativistic jet electrons Compton 
scatter an external radiation field ( ISikora et al.l Il994t 
iDermer & Schlickeiseill2002l) from the accretion disk or 
the dusty torus surrounding it. In the first case, the disk 
generates UV seed photons which are then reflected to- 
ward the jet by the broad line region within a typical 
distance from the accretion disk of the order of 1 pc. 
In the second case, the dusty torus could provide IR 
seed photons that are emitted at larger distances from the 
jet. In any of these ERC models, the VHEGR emission 
is expected to be very steady, implying long-lived TeV 
blazars, resulting in more patchy blazar heating. 

We leave detailed studies of the inhomogeneity in the 
blazar heating rate resulting from intrinsic inhomogeneity and 
variability in the spatial distributions of the blazars themselves 
for future work. 
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Figure 6. cilogQ/dlog^gL as defined in Section P^. 3. 21 for the TeV blazai's 
responsible for lialf of the total heating rate for a number of redshifts (top: 
from violet blue to light-blue: Zo = (soUd), 0.1 (long-dash), 0.3 (short-dash), 
0.5 (long-dash-dot), and 1.0 (short-dash-dot): bottom: from green to dark-red: 
Zo = 2 (solid), 3 (long-dash), 4 (short-dash), 5 (long-dash-dot), and 6 (short- 
dash-dot)). For reference, L(I)b(Zo,L-) is shown in the inset (in comoving units). 

3.3. Properties of the Blazars Responsible for the Heating 

In our discussion of the homogeneity of blazar heating we 
have necessarily attempted to define a class of objects respon- 
sible for the bulk of the heating. This is most directly done by 
considering those TeV sources which produce, say, 50% of the 
total heating rate at a given observer redshift. Given this pop- 
ulation, we may also now address the properties of the most 
relevant sources themselves; i.e., we can identify which types 
of TeV blazars dominate the heating. As we have already men- 
tioned, this cannot be too different from the set of blazars which 
dominate the luminosity density. We address this question in 
some detail in the Appendix, including constructing a simple 
analytical toy model. 

Figure|6]shows the luminosity distribution of the heating-rate 
integrand once integrated over distance and assuming our form 
for (I>b{z,L) (see the Appendix for a precise definition). The 
upper-half of the heating rate is dominated by sources with lu- 
minosities of approximately 10'*^ ergs"' at low redshifts, and 
rises to 10'"' ergs"' by z ~ 2 before declining slightly there- 
after. This evolution is driven both by the changing shape and 
the changing normalization of (ps with redshift. At all redshifts 
this is larger than the median luminosity. Thus, generally it ap- 
pears that the heating is due predominantly to high-luminosity 
objects. Due to the relatively small Dpp, these are also neces- 
sarily nearby. 

3.4. Other Potential VHEGR Sources 

Up to now, our focus on the TeV-blazars is motivated by the 
fact that the vast majority of extragalactic VHEGR sources ob- 
served are blazars. Other potential source of VHEGR emission 
exist and so we will discuss these sources in this section. In as- 
sessing their potential contribution in VHEGRs, however, we 
find that these VHEGR sources are subdominant. 
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3.4.1. Starburst Galaxies 

Starburst galaxies are characterized by the presence of rapid 
star formation, and hence, are sites of numerous supernovae. 
Supernovae are known to be efficient particle accelerators, and 
are presumed to be the primary source of the galactic cosmic 
rays. As cosmic rays can produce VHEGRs, starburst galaxies 
are potentially bright VHEGR sources. 

While many starburst galaxies contain active AGN, that 
may also emit VHEGRs, unambiguous VHEGR emission 
form the starburst itself has been seen in at least two 
cases, M82 and NGC 253, with E > lOOGeV luminosi- 
ties of 5 X 10^^ ergs"' and 2 x 1 0^^ ergs"', respectively 
([VERITAS Collaboration etani2009l: lAcero et al.l l2009). The 
luminosity of these two sources are dwarfed by the local TeV 
blazars, and thus the VHEGR flux from M82 and NGC 253 
does not contribute to the heating of the local IGM. In addi- 
tion, their VHEGR luminosity produces a pair beam which is 
too dilute for plasma processes such as the "oblique" instabil- 
ity to beat cooling via inverse Compton. Hence, they do not 
contribute to the local VHEGR flux. 

However, very high star-forming systems such as the ultra- 
luminous infrared galaxies (ULIRGs) where Ljr > lO'^L©) 
and LiR is the bolometric infrared luminosity (5 /im < A < 
1000 /im), may evade this constraint. If the VHEGR emis- 
sion from starburst galaxies is due to cosmic rays accelerated 
by supernovae, the VHEGR luminosity above lOOGeV, L is 
then proportional to the star formation rate. For starbursts, 
this is linearly r elated to the continuum infrared luminosity 
(lKennicutdll998h . Thus, normalizing by M82 and NGC 253, 
we have 

While this relationship is extremely uncertain (the normaliza- 
tion varies by a factor of two between M82 and NGC 253), 
it suggests that ULIRGs, or perhaps hyper-luminous infrared 
galaxies (HLIRGs, Ljr > IO'-'L©) may be sufficiently bright to 
contribute to the heating of the IGM. 

While the fluxes from the brightest ULIRGs remain much 
smaller than those associated with the typical TeV blazars, 
there are many more starburst galaxies than AGN. At all red- 
shifts ULIRGs constitute the high-lumin osity tail of the star- 
forming galaxy luminosity funct ions (Le Fl oc'h et al.l 120051 : 
iCaputi et al.ll2007t iMagnelU et al.ll2009: Goto et al.ll2010l) . In 
the present epoch, the density of ULIRGs is roughly 4 x 
10"^Mpc"^, and thus the corresponding VHEGR luminosity 
de nsity of these objec ts, ^ 2 x 10^''' ergs"' Mpc"-' (see Table 8 
of iCaputi et al.ll2007b . is negligible in comparison to the TeV 
blazars, roughly 5 x 10^^ ergs"' Mpc"^. 

However, due to the steep decline of the luminosity function 
in the ULIRG range, small changes in the location of the break 
luminosity result in large changes in the comoving luminosity 
density of ULIRGs. As a consequence, the comoving luminos- 
ity density of ULIRGs grows much faster than that of quasars 
(and thus presumably the TeV blazars). Nevertheless, even at 
z = 2, roughly the redshift of peak star formation, the comoving 
number densit y of ULIRGs rema ins below ~ 2 X 10"4Mpc"^ 
(see Table 8 of ICaputi et al.ll2007l) . corresponding to a comov- 
ing VHEGR luminosity density of < lO^'^ergs"' Mpc"^^, more 
than two orders of magnitude smaller than the contemporane- 
ous TeV blazar population. For this reason, we neglect the star- 
burst contribution to heating the IGM here, though they may 
represent an secondary source class. 



3.4.2. Magnetars & X-ray Binaries 

Stellar-mass objects, such as magnetars and X-ray binaries, 
generally have difficulty reaching the flux limits required for 
plasma cooling to dominate the pair beam evolution. Even 
at Eddington-limited VHEGR luminosities, reaching isotropic- 
equivalent luminosities of 10"*^ ergs"' requires beaming factors 
of roughly lO"*, corresponding to jet opening angles of roughly 
2°. In practice, the formation of radio jets is believed to be 
associated with substantially sub-Eddington accretion flows in 
these objects, and thus exacerbating the beaming requirement. 
More importantly, these sources may suffer compactness prob- 
lems: the VHEGR emission regions in stellar-mass jets must 
necessarily be very far from the central object to avoid in situ 
pair-production. Finally, were X-ray binaries and magnetars 
generally strong, persistent VHEGR emitters with sufficiently 
large fluxes, we would expect that many would have been al- 
ready detected. 

3.4.3. Gamma-ray Bursts 

Gamma-ray bursts (GRBs) are natural candidates due to their 
large luminosities and strong inferred beaming, and we con- 
sider them as a example of the general class of energetic tran- 
sient sources. Unfortunately, little is known about the VHEGR 
emission of GRBs. Presently there is a single re port of a TeV 
signal associated with a GRB (GRB 970417A, lAtkins et al] 
,2000, 2003), though due to the large distances at which they 
can be observed and comparative rarity this is not unexpected. 
However, such high-energy emission is possible in principle, 
presumably due to inverse-Comptonization of the promgt emis- 
sion and/or X-ray afterglow (see Section VIII of iPiranI [20041 
and references therein). Fermi observations of GRBs have 
shown that for many events the Band spectrum can be extended 
to - lOOGeV (Abdo et al. 2009d b.c: Ackermann et al. 20l3), 
though in at least one case a spect ral break below lOGeV 
has been observed (GRB 090926A, I Ackermann et al.l 120 lib . 
Thus it remains unclear if in practice the high-energy emis- 
sion is attenuated within the emission region. Moreover, since 
GRBs are inherently short-lived events, the luminosity lim- 
its described in Paper I, which require the VHEGR emitting 
phase to last for a plasma cooling timescale (roughly 10~yr- 
10^ yr), are not directly applicable. A similar analysis, obtained 
by limiting the beam growth time to the GRB duration, gives 
VHEGR isotropic-equivalent energies of lO^'^erg. This is com- 
parable to the total prompt and afterglow emission for only the 
bright est bursts, com prising roughly 5% of GRBs observed by 
Swift dGehrels et al.l l2009). Nevertheless, even assuming that 
all GRBs produce the requisite high energy emission, at an op- 
timistic present local rate of roughly 0.5 Gpc""* yr"' , produces a 
comoving luminosity density of < 10^*'ergs"' Mpc""*, roughly 
three orders of magnitude less than that due to TeV blazars at 
z= 1. 

4. THE THERMAL HISTORY OF THE IGM 

The previous sections have shown that the VHEGR emission 
from luminous TeV blazars heats the IGM, quantified the mag- 
nitude and stochasticity of this heating locally, and estimated 
its evolution as a function of redshift. We are now in a posi- 
tion to discuss its impact on the thermal history of the IGM 
in detail. In the following we will show that TeV blazar heat- 
ing can be substantial, dominating late time photoheating, and 
that its uniform nature naturally imprints its signature onto the 
temperature-density relation of the IGM. 
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The canonical history of the IGM is shaped by two important 
events: H reionization by the first star s at z ~ 6-10 and Hell 
reionization by quasars at z ^ 3 (see iFurlanetto & Ohll2008l 
and references therein). Hydrogen and Hell reionization both 
heated the IGM to a T - 2 x lO'^K-S x lO'^K or higher (in the 
case of Hell). Subsequently, the Universe cooled via adiabatic 
expansion, balanced by continuing photoheating due to ioniza- 
tion of recombining H. The entire canonical history of the IGM 
can thus be summarized as a competition between photoheat- 
ing and adiabatic cooling, punctuated by intervals of sudden 
photoheating 03 

The addition of TeV blazar heating adds an additional ex- 
tended heating component, and fundamentally alters the canon- 
ical picture for the thermal history of the IGM. In the following 
we explore the consequences of blazar heating usin g the one- 
zone model originally due to Hui & Gn edinI (119971 hereafter 
HG97) (see also iHui & Haim a n 200 3). We begin by introduc- 
ing this model in detail (Section l4~l1 l. describe thermodynamic 
consequences of the new heating contribution from blazars and 
relate this to high-z Lya measurements (Section l4~2] i. and close 
with a discussion of the implications for the local Lya forest 
(Section|431). 

4.1. One-Zone Model for the IGM 

The thermal evolution of a fluid element in the IGM is gov- 
erned by 



dQ 



dT ^ 2T dS T dT,Xi 

— =-2HT+ ^ - + — 

dt 3(1 + S)dt SX/ dt 3^B«bary lir 



(28) 



where H is the redshift dependent Hubble function, S is the 
mass overdensity, X, = «,/ «bary is the proper number fraction 
of species /, relative to the proper number density of baryons, 
"bary = ^BPa/mp, Pa IS the Critical mass density of the uni- 
verse, and dQ/dt is the heating and cooling rate of the gas. The 
heating and cooling of IGM gas is governed by four processes: 
adiabatic cooling/heating from Hubble expansion/gravitational 
collapse, H/He photoionization heating, H/He recombination 
cooling, Compton cooling, and heating from TeV blazars. The 
evolution of the proper number fraction of the various species 
is given by 



dXi 
dt 



- -X/Ti + ^ «baiy^j^/t^, 



ijk J 



(29) 



where the F, are the associated atomic rates, not to be confused 
with the beam instability cooling rates discussed in Section |2] 
Finally, we demand a prescription for the density evolution, 
i.e., the evolution of 6. For this, we follow HG97 and assume 
the Zel'dovich approximation: 

1+S = dct-' {Si J + D^i;ij), (30) 

where D+ is the linear growth factor ( iPeeblesI 119801) . The 
3x3 matrix, ip, is determined by initial conditions. The 
exact form of this matrix is irrelevant. What is important, 
however, is the probability distribution of the eigenvalues 
of this ma trix. For a Gaussia n random field, the solution 
is known jDoroshkevichl 119701); we use the formulation of 
iReisenegger & Miralda-Escuda(ll995l) . 

Reionization of the Universe occurs in two stages: first H 
is reionized at some large redshift by stars and later He is 

" We have ignored gravitational (shock) heating, which is progressively 
more important at densities larger than the mean density. 



reionized around a re dshift of z ~ 3 by quasars (see, e.g., 
IFurlanetto & Ohll2008h . To model this reionizing history, we 
adopt a sudden H photoionizing model (HG97): 



J(z) = 



Jo for z<z„ 

for Z > Zrf 



(31) 



where Zieion is the redshift of H reionization and Jq is the nor- 
malization of the reionizing radiation. While, this is not a 
realistic model of how H reionization occurs, its late time 
evolution (especially after z ~ 3) should be reasonably accu- 
rate. This is because the photoionizing bac kground is observed 
to be roughly constants at these redshifts (iBolton et aTl 120051: 
iBecker et afll2007l; iFaucher-Giguere et al.ll2008h and the late 
time temperature asymptotes to a single value. This "loss of 
memory" of the specific reionization history in the evolution of 
the IG M is typical of reionization calculations (iHui & HaimanI 
l2003h . 

To model the redshift dependence of H and He photoioniza- 
tion, we use the following spectral model for the radiation: 



Je{z) = J(z)(4- 



-1.6 r 1 for£<£Heii, 

<^ 0.0 for £ > EhcII & z > Zhc , (32) 

I 1.0 forE > ^Hell & Z < ZHe, 



where E is the energy of the photon, Eui, Euei, Eneii are the 
threshold energies corresponding to the ionization of HI, Hel, 
and Hell, zhc = 3.5 is the redshift of He reionization. The spec- 
tral index of -1.6 is typical of quasars and the spectral model 
of Eq uation ([32b is similar to the He reionization model stud- 
ied bv IFurlanetto & OhI ( 120081) . Our model differs from theirs 
in that we ignore the density dependent effects of He reioniza- 
tion, i.e., that dense regions are reionized first. However, since 
we are interested solely in the magnitude of the effect of blazar 
heating, adopting this simplified model is justified. 

The normalization of the photoionizing background in our 
model is fixed, determined by the H photoionizing rate: 



rHi=47r 



Em 



dE 



(33) 



where am is the photoionizing cross section of HI. We choose 
the normalization of Je to be Fhi = 5 x 10"'^, which is in- 
ferred (with significant uncertainty) fro m the mean absorp- 
tion of the Lya fore st (see for instance iBolton et al.l 120051 : 
iFaucher-Gig uere et a !1l2()08l) . 

Equations (|281 - |32] | constitute a complete model for the evo- 
lution of a fluid element in the IGM. We numerically integrate 
these equations using a prescription for the heating and cool- 
ing microphysics specified in the Appendix of HG97. In Fig- 
ure |7] we plot the evolution of the temperature for (5 = patch 
with Zieion = 19, 10, and 6. The solid lines are the purely pho- 
toionized models without the effect of additional heating. For 
each of these models, we also set the redshift of Hell reion- 
ization at ZHeii = 3.5, which results in a temperature jump at 
that redshift. We note that at late redshift, the three different 
(purely photoionized) reionization histories asymptote to a sin- 
gle temperature evolution, highlighting the "loss of memory" 
property that is generic to photoionization-dominated models 
CHui & Haiman 2003). 

Generally, the thermal history of a given patch depends upon 
the particulars of the (5-evolution of the patch. In Figure |8] we 
show the temperature of ~ 4000 realizations of an evolving 
patch as a scatter plot at a number of redshifts, ranging from 
z = 0.5 to z = 4. From this it is clear that the temperature-density 
relation is well approximated by a power law, consistent with 
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Figure 7. Thermal history of the a 5 = pat ch of the IGM for the numerical 
solution (solid Hues) of Equations j28t . i29i . and )30t using the prescription 
for the microphysics as specified in the Appendix of HG97. T he so hd c urve s 
are for sudden reionization histories for H and He (Equations )3U and (32)) 
for 

^^reion — 19, 10, and 6 and ZHell — 3.5 going from left to right. The dashed 
(dashed-dotted) lines show the evolution usin g the blazar luminos ity density, 
i.e., using the quasar luminosity density from <Hopkins et alJ2Q07l) to normal- 
ize the local heating rate in our standard (optimistic) models (see Figure fTOl and 
surrounding discussion). 




1+6 

Figure 8. Temperature-density scatter plot for 4000 realizations at ; = 3 
(black dots), z. = 2 (blue squares), z = i (green triangles), and z = 0.5 (red 
diamonds) when heating from TeV blazars is ignored. The shape of the 
temperature-density plot can be clearly be fit with a power law (with a pos- 
itive index) and steadily falls from z = 3.5, i.e., after He reionization. The 
con'espondence between the z = 4 and z = 2 points are simply an accident of 
choosing He II reionization at z = 3.5 and plotting the temperatures and densi- 
ties at z = 2 and 4. Namely the Universe adiabatically cools, but the injection 
of heat at z = 3.5 resets the temperature and it proceeds cooling from that point 
onward. 

those of HG97. Figure [8] represents a typical temperature- 
density relation, i.e., T-S relationship, that is typical of most 
reionization calculations. Low density regions in the Universe 
are cooler compared to high-density regions due to decreased 
recombination (and hence photoheating) and a more rapid ex- 
pansion (and hence greater adiabatic cooling). Missing from 
this simple picture are the e ffects of shocks a nd outflows from 
galaxies, i.e., feedback (see lDave et al.ll20Tol for instance). In 
addition, the temperature of the IGM is relatively cool (< 10"* 



Figure 9. Photoheating (black solid line) and blazar heating rates for the stan- 
dard (dashed line) and optimistic (dot-dashed line) as a function of redshift. 
The fits (eq.(38)) are denoted by gray lines for the standard and optimistic 
blazar heating models. We use a 5 = fluid element in the IGM and calculate 
the amount of heat per baryon (in units of eV) per Gyr The sharp jump at 
z = 3.5 is due to Hell reionization. Aside from this point, it is clear that blazar 
heating dominates the heating of the IGM. At late times z < 2, it is larger by 
over a factor of 10, though this is substantially greater for underdense regions. 

K) for S <Q. This is because the temperature is suddenly raised 
to a few X 10"* K after H and He reionization, but rapidly cools 
due to the effects of adiabatic expansion. These generic char- 
acteristics are typical of most reionization mo dels and are ex- 
pected following the arguments of HG97 and iHui & HaimanI 
il2003i) . 

4.2. Contribution of TeV Blazar Heating 

When the effects of heating due to the VHEGR emission 
from blazars are included, the properties of the IGM are sub- 
stantially altered. The consequence of TeV blazar emission is 
shown for a S = patch of the IGM by the dashed line in Fig- 
ure|2l The thermal history begins to deviate significantly from 
that due to photonionization and adiabatic expansion alone by 
Z — 6, becomes dominant near z — 3, and peaks at roughly 4- 
8) X lO'^K at z ~ 1 before the rapid decline in Ab(z) combined 
with adiabatic cooling causes the temperature to fall off (the 
range corresponds to the uncertainty in estimating the number 
of blazars contributing to the heating rate). Thus, it is clear 
that heating by blazars is significant, dominates at low red- 
shifts (following He reionization), and potentially dominates 
the thermal evolution of the IGM in low-density regions. 

The effect of TeV blazar heating qualitatively changes the 
picture of the IGM. First, the temperature-density relation is 
inverted with the low-density regions being the hottest. Sec- 
ond, the overall temperature of the IGM is significantly hotter. 
The reasons for both of these are twofold. First, TeV blazars 
are a substantial reservoir of heating, potentially increasing the 
IGM temperature by a few x lO'^K, and dominating the con- 
tribution from ionizing photons for l+S < 10. Second, the 
heating rate is nearly independent of density, depending most 
strongly upon the number density of TeV blazars and the num- 
ber density of UV photons, both of which are nearly uniform 
(though see Section [32]). The effect of a uniform heating rate 
is that the energy deposited per baryon is substantially larger in 
more tenuous regions of the Universe, with underdense regions 
experiencing larger temperature increases as a result. 
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The dominance of TeV blazar heating over photoheating is 
shown explicitly in Figure |9] where the blazar heating rate 
(dashed line) is compared to the photoheating rate (solid lines) 
as a function of redshift for a (5 = patch of the IGM. At z = 3 .5, 
there is a sudden jump in the photoheating rate due to nearly 
instantaneous Hell reionization. Following Hell reionization, 
the blazar heating rate is about an order of magnitude larger 
than that due to photoheatin^3 Because the photoheating rate 
is oc (1 + (5), the dominance of TeV blazar heating is even more 
apparent at lower densities. 

The contribution of TeV blazar heating to the thermodynam- 
ics of the IGM for the standard model appears to be signifi- 
cant around the period of Hell reionization. In our model, this 
is partially a result of our sudden reionization prescription for 
Hell reionization at z = 3.5. However, we can also show with 
the following order-of-magnitude calculation that the effect of 
TeV blazar heating must begin to be important around the era 
of He II reionization that has been observationally constrained 
to be around z ~ 3. 

To begin we first show that He reionization finishes around 
z ~ 3. The comoving number density of He is 



«He = fl 



He 



Amnip 



1.5 X lO'^^cm" 



(34) 



where /hc = 0.24 and Ane = 4 are the primordial mass fraction 
and the atomic number of He, respectively. To estimate the 
comoving density of Hell ionizing photons at z ^ 3, we note 
that t he I Ry photon comoving density at z = 3.5 from quasars 
in the lHopkins et al.' (2007) QLF is /iiRy w 5 x 10"^cm"^Gyr"' 
(e.g. see Figure 9 of Hopkins et al. 2007). As the spectral index 
of quasars is -1.6, this implies that comoving number density 
of ionizing Hell photons is ri4Ry w 5 x 10~**cm"^Gyr"'. Thus, 
the total comoving density of Hell ionizing photons produced 
at z = 3.5 is 



n4Ry ' 



Hiz= 3.5) 



10 



-7 -1 

cm . 



(35) 



Before we compare of equation dJSl ) with ( l34b . we note the it 
takes roughly 2-3 He ionizing photons to completely ionize He 
and that obscurin g material around a qu asar will remove half of 
the ionizing flux jMcOuinn et al.ll2b09l) . With these efficiency 
factors in mind, the comoving density of Hell ionizing photons 
at z ~ 3 is just large enough to reionize Hell. 

The amount of excess energy per Hell ionization is 16 eV for 
an ionizing radiation spectral index of - 1 .6. Given that there are 
w 3 ionizing photons per Hell reionization, the excess energy 
dumped into the IGM during Hell reionization is 



GexcHell - fexcHeU 



CI 



:6x 10"^eVcm"^ (36) 



where SexcHeii = 48 eV is the average excess energy dumped 
per Hell ionization, fm = 0.24 is the primordial mass fraction 
of He, and Ahcii = 4 is the atomic number of He. To compare 
this to the energy dumped from TeV blazars, we start with the 
comoving luminosity density of blazars at z = 3.5 is Ab{z = 
3.5) « 3 X lO-'^ergss"' Mpc"-', which we estimate from Figure 
[3] Over a Hubble time at z = 3.5, the amount of energy dumped 

Prior to Hell reionization, tlie pliotoionization and blazar lieating rates are 
inco nsis tent due to the artificial ionizing photon distribution assumed in Equa- 
tion 1321 Specifically, for reasons of simplicity, we have ignored the ionizing 
photons produced by the quasars. Following Hell reionization, however, this 
is no longer an issue. 



into the IGM by TeV blazars per comoving volume is then 
Ab(z = 3.5) 



Q 



B|;=3.5 = 
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■ 5x 10"^eVcm"^ 



(37) 



which compares favorably with Equation ( l36l ) and shows that 
as Hell reionization is being completed, the effect of blazar 
heating begins to be pronounced. 

For convenience, we provide a fitted third-order polynomial 
formula for our estimated proper blazar heating rate for z < 
5.7: 



log 10 



leVGyr- 



= 0.0315(l+z) -0.512(l+z)'' 



+2.27(1 +z)-logioe„od 



(38) 



Here, logigQmod = {2.38,2.08} for the "standard" and "opti- 
mistic" (see below) blazar heating model, respectively. The fits 
are shown by the solid (gray) lines in Figure |9l Note that we 
have calculated the heating rate as a heating rate per baryon, 
i.e., 2/«bary, SO that it is independent of the ionization fraction. 

The effect of blazar heating upon the void temperature- 
density relation is shown in the left panels of Figure[TO]for the 
"standard" (upper) and "optimistic" (lower) heating models. In 
comparison to the case in which blazar heating is neglected 
(Fig. |8]l, the voids are hotter at all redshifts. Overall, the IGM 
never drops below lO'^K, with low-density regions remaining 
substantially hotter than this. By contrast, the temperature of 
the IGM without blazar heating does not exceed 10"* K (Fig. [8|. 
In these low-density regions the temperature can exceed 10^ K 
with TeV blazar heating, almost two orders of magnitude hotter 
than anticipated by photoheating and adiabatic cooling alone. 
To show this result more clearly, we show in the right panels 
of Figure[TO]the ratio between the temperature for a photoheat- 
ing only model (Tphoto) and for a model including blazar heat- 
ing (Tbiazar) for the "standard" (upper) and "optimistic" (lower) 
heating models. This shows that the temperature in low-density 
patches in the Universe, i.e., small l + S, increase by a factor of 
^ 100 and ^ 200 for the "standard" and "optimistic" heating 
models, respectively. 

Figure [TO] shows that blazar heating qualitatively changes 
the IGM temperature-density relation. In particular, the IGM 
temperature-density relation is inverted for S <0. This com- 
pares favorably with recent evide nce for jus t such a n inverted 
tempe rature-density relation (B olton et al.l 120081 : IViel et all 
120091) . Notably, in the lower panel of Figure [lO] we plot as 
a soUd (gray) line the temperature-density relation T = 2.4 x 
10^(1 + S)-°-^'^K, inferred at z = 3 empirically bv IViel et al.1 
(120091) . Here to produc e a bet ter match to the empirical mea- 
surements of IViel et al.l (l2009h . we introduce the "optimistic" 
heating model, i.e., we choose a value of the remaining co- 
efficient of the sky incompleteness correction for TeV blazars 
of 7]sys = 1 -6, which includes variations of the TeV blazar red- 
shift evolution, additional TeV sources that may contribute to 
the plasma instability heating, and potentially spectral variabil- 
ity. As a result, the local heating rate is enhanced over the 
"standard" heating rate by a factor of two or ~ 1-4 x 

10~^eV cm ~^Gyr~' . The mat ch produced between the empirical 
relation of IViel et al.l (l2009h and our "optimistic" TeV blazar 
model is striking in that it matches the slope in the absence of 
any tuning as well as the normalization within the uncertain- 
ties of our estimated incompleteness correction (see Figure[3]). 
This inversion of the temperature-density relation is difficult to 
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Figure 10. Temperature-density scatter plots (left panels; same as Figure [8) including the effects of blazar heating for our standard (top) and optimistic (bottom) 
heating models. The right panels shows the ratio between including only the effects of photoheating (Tphoto) and including the effects of blazar heating (7i,iiiz„) for 
the same realization for the standard (top) and optimistic (bottom) blazar heating models. Note for the left panels that the temperature evolution due to blazar heating 
is qualitatively different from the earlier case (Fig. [8). Starting at ; = 4, the temperature-density relation in the lowest density regions of the Universe is inverted. This 
inverted temperature relation moves to higher densities with decreasing redshift until aX z = 0.5, it goes down to 5 2. The IGM temperature is also significantly 
hotter , with the hottest regions having T > 10^ K. Overplotted is the observational determined temperature-density relation at z = 3 (solid gray line) from jViel et si\ 
120091) . The right panels show that the temperature increases by a factor of ~ 100 (standard heating) and ~ 200 (optimistic heating) for small 1 + 5 patches of the 
Universe when the effects of blazar heating are taken into account. 



reproduc e using Hell reionization alone (iMcOuinn et al.ll2009l: 
iBolton & Becker 2009) though is a natural consequence of TeV 
blazar emission. 

Energetically, the magnitude of the impact from TeV blazars 
in low-density regions is somewhat surprising. The radiative 
output from stars and quasars far outstrips that from VHEGR 
sources, yet in practice the heating rate from blazars is much 
larger. This is because the photoheating rate is ultimately lim- 
ited by the HII recombination rate. Indeed, this is precisely the 
property invoked to show that the effect of the epochs of H and 
Hell reionization is washed out at low redshift, unless these 
occurred recently (iHui & Haimanll2003l) . 

By contrast no such limitation exists for TeV blazars, paving 
the way for them to dominate the heating of the low-redshift 
Universe, a point that we have made already in the Introduc- 
tion. Instead, the VHEGR emission from blazars is deposited 
in the IGM with order unity efficiency. The rate of heating in 
this case depends linearly on the radiative output of these TeV 
blazars, and is independent of the atomic physics of the IGM. 
In addition these sources have a constant heating rate per unit 



volume (whereas photoheating heats per unit mass). Hence, the 
effect of blazar heating, which is already pronounced, is ampli- 
fied relative to photoheating in low-density regions, leading to 
an inverted temperature-density relation. Therefore, not only is 
the memory of photoheating is erased in low-density regions, 
but it is overwritten by the record of blazar heating. 

Following this work, we have recently completed a study 
of the effects of blazar heating in a more detailed hydrody- 
namic model of structure fo rmation, reported in a follow-up 
paper (iPuchwein et al.ll201 lb . Using the blazar heating pre- 
scription given by equation dSST l, we show that the optical depth 
of the Lya forest is reproduced using a H photoionization rate 
of Phi ~ 5 X 10~'-'s"' or equivalently using the inferred ioniz- 
ing background from (jpaucher-Giguere et al. 2009). We also 
confirm that the low density IGM again possesses an inverted 
temperature-density relation (as shown in Figure [Tol i. In ad- 
dition, a detailed comparison between the results of our nu- 
merical calculations and observations show that a blazar heated 
universe matches the one- and two- points statistics of the high 
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redshift Lya forest, but also matches the line width distribu- 
tion. This excellent agreement was achieved using the best es- 
timate on the evolution of the photoionizing background (with- 
out tuning) and is due to the excess energy injection into the 
low density IGM. 

4.3. Implications for the local Lya forest 

While TeV blazar heating can with significant success repro- 
duce some of the peculiar properties of recent high-z Lya ob- 
servations, the most dramatic departures from the IGM thermal 
history in the absence of blazar heating occur at z < 1 . This is 
a result of the relatively recent nature of quasar activity and the 
cumulative effect of blazar heating. Thus we might anticipate 
dramatic consequences for the local Lya forest, potentially in 
conflic t with studies of nearby Lya c louds (e.g.. lPenton et al.l 
l2000bllaLl2002ll2004l:lDave et al.l2010l) . However, this does not 
occur due to the physical nature of the clouds that produce the 
local Lya forest. 

The structures responsible for the local Lya forest are almost 
certainly associated with significant overdensities a s suggested 
by l a rge -scale hydrodynamic computations (e.g., iDave et al] 
120 lO r^. These simulated Lya forest calculations provide a 
means to relate the empirically inferred HI columns to the 
properties of the dynamical structures responsible for the Lya 
clouds. For HI columns of 10'^-10'"*cm"^, the lowest values 
for which Lya measurements exist, the simulations find that 
the local Lya forest is produced primarily by the intergalactic 
filaments, corresponding to 1 +(5 > 5- 500 (e.g., see the b ottom- 
left panel of Figure 8 and Figure 9 of iDave et ani2010l) . 

We emphasize this point with the following order-of- 
magnitude estimate. An upper limit upon cloud sizes can be 
obtained using the line width s, b, translating int o a proper 
size of L4(/7/10-kms"')Mpc (iPenton et al.ll2000al) . A lower 
limit upon the local ionizing background can be inferred from 
the optical depth of VH EGRs, giving ~ 10~^ergs"' cm"^sr"' 
(lAharonian et al.l l2006ah . Employing ionization balance to 
set the HI fraction, for z < 1 this gives the column densities 
<3 X lO'^l +(5)2(1 +z)^-'^(r/10'*Kr"''5 forr < lO^K. From 
this it is clear that the nearby Lya clouds must correspond to 
regions with 1 + 5 ^ Ipl 

At these overdensities, the impact of blazar heating on their 
thermal history is modest at best. By 1 + 5 = 5 the blazars 
change the IGM temperature by a factor of 2, raising it to 
roughly 4 x 10"* K, co mparable to the t emperatures typical of 
flie local Lya clouds (iDave et al.l 1201 Oh . At 1 +i5 > 50 blazar 
heating is negligible in comparison to photoionization and 
shock-heating. Thus, because low-z Lya absorbers are bi- 
ased toward high-density regions, they serve a poor probes of 
low-density regions and are unaffected by the effects of TeV 
blazars. 

The inclusion of blazar heating into these large-scale hydro- 
dynamic calculations is an important next step. However, we 
do not expect a significant impact on the local Lya forest. A 
more important impact of TeV blazar heating would be on the 
formation of collapsed structures, which is a topic that we will 
explore in Paper III. 

4.4. Limits on the direct detection of the hot IGM 

We note in passing, that for z > 2, the same simulations imply that Lyo 
forest contains significant contributions from regions with 1 + 5 ~ 1 . 

Note, however, that since the column density is oc (1+;)'' atz > 1, we 
would nevertheless expect low-density regions to contiibute significantly to the 
high-4 Lya forest, as has indeed been found to be the case. 



As low-density regions are immune to local Lya probes, we 
now turn to the question whether there are methods of directly 
inferring the presence of such high temperatures in these re- 
gions. First, we compute the mean Comptonization of the 
CMB due to blazar heating, 

{y)=aT I dD 5 , (39) 

Jo nieC- 

where iie is the free electron density. Here, lie is the physical 
density of free electrons and we integrate along the proper dis- 
tance of the photons back to recombination, dD = -cdad~^ = 
-cda{aH)~^ where the minus sign arises due to the choice of 
coordinates which are centered on the observer Performing the 
integral with our temperature evolutions at mean density (i.e. 
neglecting gravitational heating by formation shocks) yields 
mean Comptonization values of (y) = {1.4, 1.9,2.5} x 10"^ for 
our models without blazar heating and those with standard and 
optimistic blazar heating. To date the best limits on the mean 
Comptonization come from the COsmic Background Explorer 
Far-InfraRed Absolute Spectrophotometer experiment (COBE 
FIRAS) which measure the dif ference between th e CMB and 
a perfect black-body spectrum ( iFixsen et al.lll996l) . Their up- 
per limit of \y\ < 1 .5 x 10"^ (95% confidence level) is perfectly 
consistent with our inferred mean Comptonizations. With cur- 
rent technology, it appears to be quite feasible to measure the 
deviation of the CMB spectrum from a perfect blackbody form 
with an accuracy and precision of 1 ppm yielding constraints 
on the cosmic y-parameter at the lev el of 10"^ and provide a 
spectrum of the anisotropy to 10% (IFixsen &Matheril2002h . 
We note, however, that the mean Comptonization is expected 
to be dominated by gravitational heating which leads to values 
of (y) = 2.6 X 10"^ in ferred from cosmological simulations by 
ISpringel et al.l (120011) . The mean temperature of these cosmo- 
logical simulations of (To) = O.BkeV suggests that the signal 
is dominated by collapsed galaxy groups and unvirialized in- 
fall regions onto galaxy clusters that constitute the hot compo- 
nent of the warm-hot intergalactic medium. Hence, in order to 
measure the blazar heating signal in the mean Comptonization 
of the CMB, the fluctuating part due to gravitational heating 
would have to be subtracted first. Since galaxy groups domi- 
nate the fluctuation power on angular scales of < 5', the mea- 
surement of the deviation of the CMB spectrum would have to 
be performed on these small angular scales which appears to 
be difficult. 

Second, we estimate the emission of the IGM due to free-free 
emission (bremsstrahlung) and synchrotron emission. Turn- 
ing to the question of synchrotron emission, we estimate the 
synchrotron frequency using a maximal IGMF strength of B 
10"^ G. At temperatures of ^ 10'*- 10'' K, the electrons in the 
IGM are nonrelativistic. Hence the synchrotron frequency is 
roughly the Larmor frequency, Wsync ~ ojg = eB/nieC w 2 x 
10~2 (B/10~^G) s"'. These low frequencies are well below 

the plasma frequency of the IGM, LOp = ^/ATieh^e/nie ~ 25 Hz 
and, hence, will be absorbed by the IGM. Free-free emission 
from the IGM c an cause a distortion in the CMB at low fre- 
quencies (iBartle tt & Stebbins 1991). As the free-free emis- 
sivity scales like oc HgT""^^ for fully ionized primordial gas, 
it scales weakly with temperature but strongly with dumpi- 
ness. As a result, overdense haloes dominate free-free emis- 
sion, with the contribution fr om low density IGM being 1-3 
order s of magnitude smaU er (lOhl 119991: ICoorav & Furlanettol 
12004 lPonenteetani2011l) . The effect of blazar heating will 
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further suppress this free-free emission (albeit mildly) from the 
IGM compared to that of overdense halos. Finally, the best 
constraints on this fre e-free optical depth of Tff < 1.9 x 10"^ 
jBersanelli et al.l ll994') is still too large compared to optimistic 
models of free-free emission from overdense haloes by at least 
an order of m agnitude and by R:i 3 orders of magnitude from th e 
smooth IGM (ICooray & Furlanettol2004tlPonente et alJl201 Ih . 

5. MEASURING THE HIGH-ENERGY LUMINOSITY OF BLAZARS 

Thus far, we have concerned ourselves with the impact of 
TeV blazars on the Universe at large. In this section, we in- 
vert this argument and discuss how the Universe at large can 
be used as a probe of the physics of VHEGR photons and the 
global properties of TeV blazars. Namely, the estimated heat- 
ing rates used in this paper - standard and optimistic - suffer 
from a number of uncertainties, including statistical fluctua- 
tions in the number of presently observed VHEGR sources, the 
redshift evolution of the TeV blazar luminosity density, and 
the detailed form of the VHEGR and EEL spectra. However, 
due to the high efficiency with which the VHEGR luminosity 
of blazars is converted to heat within cosmic voids, the tem- 
perature history of the voids themselves provides a way to de- 
termine the cumulative TeV blazar luminosity density empiri- 
cally. 

This is possible due to a number of fortuitous properties of 
voids and TeV blazar heating. First, the high efficiency with 
which the VHEGR emission of blazars is converted to heat 
implies that even in voids the Universe acts as a calorimeter 
Second, blazar heating dominates the thermal history of voids 
for z < 4 (see, e.g.. Figure |9]l so that this calorimeter is un- 
contaminated. We caution that Hell reionization complicates 
matters somewhat, though in our model makes a comparatively 
small correction to the temperature evolution of low-density re- 
gions (see Figure|7ll. However, this depends upon the particu- 
lar manner in which Hell reionization occurred, and could in 
principle provide a somewhat larger contribut ion in inhomoge- 
neous reionization scenarios (see Section 3 of iFurlanetto & OhI 
I2008L for a detailed discussion). Finally, within voids the adi- 
abatic losses to Hubble expansion are very accurately modeled 
in the linear regime, substantially simplifying the interpretation 
of void temperature histories. 

Perhaps the greatest uncertainty is the physics of the mech- 
anism which converts pair beams from TeV blazars to thermal 
energy in the IGM. In the context of the "oblique" instabil- 
ity, what is unclear is the luminosity cutoff of TeV blazars 
probed in this manner, which arises from the competition be- 
tween inverse-Compton and plasma processes. If this cutoff 
may be conservatively assumed to lie at an isotropic -equivalent 
ELe — lO'^^ergs"', corresponding to a true luminosity 2-3 or- 
ders of magnitude lower due to the presumed jet beaming, then 
this is dimmer than all but two of the TeV blazars known. As- 
suming that the TeV blazar luminosity function follows that of 
quasars, the luminosity density is dominated by sources near 
the break luminosity, and thus as long as the lower-luminosity 
cutoff is sufficiently low (below that of the break luminosity) 
it may be neglected. For the TeV blazars listed in Table [1] this 
break occurs near 3 x 10"*'*ergs~', and thus is well above the 
relevant cutoff. Measurements of the thermal history of voids 
directly corresponds to the TeV blazar bolometric luminosity 
evolution. In this manner, the thermal histo ry of cosmic v oids 
provides an analogous argument to that by ISoltanI (119821) for 
determining the blazar luminosity density, Ab(z). Such a con- 
straint on the history of Ab(z) can be used to study the history 



of accretion in the Universe and the jet forming efficiency by 
comparing it to the history quasar or active galaxy luminosity 
density. While probes of the low-density IGM at low redshift 
are sparse, the situation at z ^ 2 -4 is much more hopeful as we 
have discussed in Section l4~2l Thus, precision measurements 
of the Lya forest at z ^ 2-4 (e.g.lViel et al. 2009) al ongsid e de- 
tailed studies of Hell reionization ( iMcOuinn et afl 120091 e.g.) 
offer the best constraint both on the fate of VHEGR photons 
and the evolution of TeV blazars, i.e., the blazar luminosity 
density. 

6. CONCLUSIONS 

In this work, we have explored the effect of TeV blazar heat- 
ing on the thermal history of the IGM. We have argued that 
VHEGRs that are sufficiently hard to pair produce off of the 
EEL, will inevitably dump the majority of their energy into the 
IGM via plasma instabilities. Ey collating the nearby 28 TeV 
blazars with firm spectral measurements, we have determined 
the local observed heating rate after correcting for the various 
selection effects using Fermi observations of TeV blazars. This 
local heating rate of g = 7 x 10~^eVcm"^Gyr"', which we call 
the standard mod el, can be extended to higher redshift by nor- 
malizing it to the lHopkins et al.l ( 120071) quasar luminosity den- 
sity. This follows from the important result of Paper I, which 
shows that the local observed blazar luminosity function is well 
in line with the local quasar luminosity function corrected by 
a factor of w 10"^. This TeV blazar heating should be rela- 
tively homogeneous at all redshifts z < 4 with greater spatial 
variations for at higher z, approaching order unity at z ~ 6. 

This redshift dependent blazar heating is substantial and is 
larger than the photoheating rate by a factor of 15 (standard) 
- 30 (optimistic) after He reionization for a (5 = patch of 
the Universe. Using a simple one-zone model of the IGM, 
we demonstrate that the effect of including TeV blazar heat- 
ing versus not including TeV blazar heating leads to quali- 
tative and quantitative changes in the thermal history of the 
IGM. First, the injection of heat into the IGM by blazars sub- 
stantially increases the temperature of the IGM. In the case 
with blazar heating, the temperature of the IGM stays above 
10"* K, with some regions approaching 10^ K, whereas with- 
out blazar heating, the temperature of the IGM tends to stay 
below 10"* KH Second, the even volumetric heating rate of 
blazars impacts the thermal history of low-density regions 
much more strongly than higher density regions. Low-density 
regions are substantially hotter as a result with temperatures 
in excess of > 10^ K. Higher density regions on the other 
hand are not heated as much by blazars. This naturally pro- 
duces an inverted tempe rature-dens i ty rela tion which matches 
the empirical results of iViel et alj ( 120091) if we increase the 
amount of blazar heating, i.e., the optimistic model, io Q = 
1.4 X 10~^eVcm"^Gyr"'. It also provides an encouraging en- 
dorsement of our model as such inverted temperature-density 
relations are difficult to produce in standard reionization histo- 
ries. We have demonstrated these salient points more explicitly 
in a follow up paper (Puchwein et al. 2011), where we calcu- 
late the effect of blazar heating in a hydrodynamic realization 
of the universe. We show that the comparison between a blazar 
heated universe and observation of the high redshift Lya forest 
gives excellent quantitative agreement. 

As our model predicts a substantially hotter low-density IGM 

" Again, we have ignored the effects of gravitational (shock) heating which 
is important at densities above mean density. 
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that standard models predict, we then proceeded to investigate 
either if this model breaks current constraints on the local tem- 
perature of the IGM or can be directly measured. Unfortu- 
nately, the local Lya forest is an ineffective probe of this envi- 
ronment compared to the high-z Lya forest as the regions that 
give rise to the local Lya forest are dense regions that remain 
relatively unaffected by the effects of blazar heating. Other 
means of directly probing this hot IGM via Comptonization of 
the CMB and free-free emission emission are also unlikely. 

Finally, we note that the thermodynamics of the IGM can 
be used as a calorimeter for VHEGR emission of blazars in 
the Universe. Namely, because the low-density IGM is so sen- 
sitive to the total amount of energy dumped into it, which is 
dominated by TeV blazars, we argue that the thermal history 
of the low-density IGM can be used to measure the total en- 
ergy output in VHEGRs over cosmic time. In principle, this 
would allow a determination of the blazar luminosity density 
as a function of redshift, as well as constrain the history and 
physics of accretion onto supermassive black holes, i.e., rates 
of radiative versus radiative inefficient accretion and jet forma- 
tion efficiency. However, the contaminating effects of He II 
reionization would have to be explored in detail before such 
physics can be elucidated. 
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APPENDIX 

A. DEFINING THE NUMBER OF BRIGHT BLAZARS 

Here we provide a more detailed discussion of how "likely" 
it is for a patch of the Universe at a given redshift to expe- 
rience a large deviation from the average TeV blazar heating 
rate. We will pursue this primarily by attempting to compute 
the "number of blazars a patch sees" as a function of redshift, 
A/b(z). This is, however, a poorly defined quantity, the primary 
difficulty being the determination of which objects to count. 
In principle, we would like to take a census of those sources 
"responsible for the bulk of the heating"; in practice this is an 
ambiguous definition. Thus, here we will describe and contrast 
a number of possible definitions. In the process, we will also 
elucidate which objects dominate the heating. 

A.l. Preliminary Definitions and Assumptions 

In the interest of completeness we will forgo the assumption 
that the heating is local, and perform the appropriate cosmo- 
logical calculations. However, since we ultimately would like 



100 




l + z 



Figure 11. Same as Figure \5\ with the addition of the cyan, short-dash- 
dotted line, showing the instantaneous number of blazars within a comov- 
ing volume of radius (1 +z)Dpp, and the green, long-dash-dotted line, show- 
ing the number of blazars seen at a particular redshift above a flux limit of 
4 X 10"'^ ergcm"^ s"' , roughly that inferred from the very-high energy gamma- 
ray observations. The remaining lines are as defined in Figure [3] The blue, 
short-dashed line shows the instantaneous number of blazars with intrinsic 
isotropic-equivalent luminosities above 10*^ ergs"' within the r = 1 surface, 
the red, dotted line shows the number of blazars above the median luminosity 
at each redshift within the f = 1 surface, and the black, solid (long-dashed) 
line shows the number of blazars with individual heating rates that exceed 
that above which half (0.75 times) of the heating is produced. For reference, 
our estimate of the sky-completeness corrected number of TeV blazars ob- 
sen'ed in the TeV is shown by the black circle, with the eiTor bars denoting 
the Poisson uncertainty only. Note that for all calculations we have assumed 
Lm = 2 X lO''* ergs"' and adopted a = 3. 

to identify Mb as a function of observer redshift, Zo, we require 
observed-dependent definitions of the standard compliment of 
cosmological distances. Specifically, we use the following gen- 
eralizations of the standard proper, angular diameter, and lumi- 
nosity distances. 



Dp(z\Zo) = Dp(z)-Dp{zo) 
Dc{z)-DdZo) 



Da{z;Zo)-- 



Dl(z\Zo) = {1+z) 



l+z 

Dc(z)-Dc{z„) 



(Al) 



where Dciz) and Dp{z) correspond to the z = comoving and 
proper distances, respectively. Note that these reduce to their 
expected expressions when Zo = 0. 

While we have already defined the mean free path of high- 
energy gamma rays in Equation (O, we also require a definition 
of the optical depth that a gamma ray that originates at z with 
observed energy E at Zo accrues during its propagation: 



t{E,z;z„) = 



1 



dD[ 



Dpp[{l+z)E/{l+Zo),z] dz 



dz- 



(A2) 



As with t(E,z), this differs from the definition of te(E,z) at 
Zo = given in Paper I, where there we set E to the emitted 
energy of the gamma ray. 

We define the flux, F(E,„,Em), to be that integrated between 
a given energy range, E,„ to Em (usually 100 GeV to 10 TeV), 



Em 



F(E,„,Em)= I dEFE=foE^ I E'-'^dE, (A3) 



where the observed photon number flux is given by Equation 
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iT2\ . For a source at redshift z, this corresponds to an absorp- 
tion corrected flux of 



F{z;z„,E„„Em) = 



Em 



dEpEC 



-T(E.z;z„) 



(A4) 



and therefore an intrinsic luminosity between energies 
E„il+z)/il+Zo) and Em(1 + z)/(l + Zo) of L(z;zoX„E'm) = 
4'!tDl(z;Zo)^F(z;Zo,E,„,Em) where and E'^^ are the energies 
bounding the relevant range at z- Since we would like to com- 
pare luminosities within a band across redshifts (i.e., keep E,„ 
and Em fixed), we must correct for the fixed spectral shift in- 
duced by the different redshifts. Assuming that the spectrum is 
a power law (as we shall do in all cases here), this implies an 
additional redshift factor: 



L(z;zo,E„,Em) = 4tt 



■■4-TT 



1+z 

l+Zo 

1+z 



a-2 



De(z;Zo)F(z;Zo,E,„,Em) 



a-2 



Di(z;zo)e-^F(E,„,EM), 

(A5) 



where r is the spectrally average optical depth: 



t(z;zo) ■■ 



-In 



Em j^i-ct -r(E,: 



~"'>dE 



J^"E^-°'dE 



(A6) 



Note that this is closely related to Dpp(z„). 

Finally, we will adopt the form of the physical density of 
TeV blazars, (I)b{L,z), described in the main text assuming a 
cutoff at Lm = 2 X 10"** ergs"', consistent with physical models 
of hig h-energy gamma ray blazars (see, e.g., Ghisellini et all 
l2009h Associated with this we have a generalized number den- 
sity, 

flog 10 ^-M 



$b(z;L„,Lm) = 
and luminosity density. 



Mz,L)d\og^oL (A7) 



log|oi,„ 



login 



LMz,L)dlogioL, (A8) 



logioim 



defined within a given intrinsic luminosity range. 

A.2. Possible Definitions of Mb 

We now present various ways to define the number of high- 
energy gamma-ray blazars that are relevant for heating. These 
are compared in Figure [TT] which supplements Figure |5] with 
additional approximations for A/^. 

A.2.1. Local number within a mean free path 

Our first definition is the simplest one can imagine; choose 
an intrinsic luminosity range and use the local density to de- 
termine the number within a volume defined by the spectrally- 
average mean free path. That is. 



AlT - 

Nb.i(Zo\L„,) = yDIp(Zo)^b(Zo;L,„,Lm) . 



(A9) 



This is an approximation of the number of sources with intrin- 
sic luminosities in the specified range within an approximation 
of a single optical depth. This is shown for L„, = lO'^^ergs"' 
by the cyan dash-dotted line in Figure [TTI Typically, where the 



blazar population is rapidly evolving it tends to be a poor esti- 
mate of the number of sources, overestimating this number by 
nearly an order of magnitude at z > 3. 

A.2.2. Number within r = 1 

As long as (f>B{z,L) evolves slowly and Dpp <C c///o it is un- 
necessary to perform the relevant redshift integral to get the 
volume element. However, this is not always the case, and 
thus a more accurate estimate is obtained by integrating the 
blazar number density within the volume specified by unit op- 
tical depth. That is. 



NB.n{Zo\L,„)-- 



' AnDl{z;zo)^^B(z;L,n,LM)dz, (AlO) 
az 



where z\ is defined implicitly by f{zi\Zo) = 1- 

For this definition, we may choose L,„ in a variety of ways. 
In Figure [TT] we show two in particular: that from setting 
L,„ = 10"*^ ergs"' (blue short-dashed line), which may be di- 
rectly compared with the case shown for Nbj, and setting 
Em = Lo,5 as defined by Equation ( |26] ), i.e., the luminosity- 
weighted median luminosity, above which sources produce half 
of the total local luminosity density (red dotted line). For 
z < 3 this approximation for Mb gives similar results to Afs.i 
for a fixed L„,. At high redshifts, where (f>B is rapidly decreas- 
ing, the two diverge substantially. The difference between the 
A/'b,//(z;Lo.5) and the previous two is more striking, and a con- 
sequence for the assumed evolving luminosity distribution of 
blazars, which is clearly important. 

A.2.3. Number of objects above a flux limit 

While a fixed intrinsic luminosity limit may be useful con- 
ceptually, even idealized surveys do not directly probe such a 
population. Instead, most surveys are flux limited, and thus we 
also define a flux-limited definition of Mb- In this case we set 
the minimum luminosity via 



im(z;Zo,/v,i) = 47r 



1+z 



Dt{z;Zo)e 



where F„, is a fixed flux limit, from which we obtain 



(All) 



J^B.!n(Zo;F„,)= I 4TrDl(z;Zof 



dDp 

j^\Z\Zo)~. 

dz 

X $B [z;L,„(z; ]dz, (A12) 

This corresponds to the number of objects a flux-limited sur- 
vey (in the lOOGeV-lOTeV band) performed by an observer 
at redshift Zo would detect. As such, it is the most directly com- 
parable to the number of TeV blazars that have been observed. 
This is shown for flux limit comparable to that inferred for the 
TeV sample in Table[Tl 4.19 x 10"'^ergcm"^s"' in Figure [TTI 
by the green long-dash-dotted line. Note that the number of 
objects found at z = corresponds nicely to the number ob- 
served after correcting for incompleteness of present TeV sur- 
veys. Of course, this is by construction since (f>B was obtained 
from the observed population. Nevertheless, it is striking that 
many more objects would have been observed above this flux 
limit during earlier epochs, vastly exceeding any of the preced- 
ing approximations of Mb- However, it does not necessarily 
follow that all of these sources will have contributed substan- 
tially to the heating rate. 
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A.2.4. Number of objects above a fractional heating rate imposed 
heating rate limit 

The most relevant approximation for Ag, and the one we 
adopt as our primary definition, is set by the heating rates di- 
rectly. The general idea is to do what we do naturally at Earth: 
arrange all of the sources visible by an observer at Zo by the lo- 
cal heating rate they induce, from largest to smallest, and count 
until a fixed fraction of the total heating rate is reached. That 
is, set Mb to be the minimum number of sources (on average) 
required to produce a given fraction of the total heating rate. 
Given its direct connection with the heating rate, this provides 
the most natural definition of the number of sources "responsi- 
ble for the bulk of the heating." To do this, however, we first 
must explicitly define the heating rate in terms of the appropri- 
ate functions. 

Given a fixed spectrum, there is a linear relationship between 
the local flux and the heating rate, defined by Equation (|9|i: 



1 = 



Em 



Fec 



-t(E.z;Zc) 



E,„ Dpp[E{l+z)/(l+Zo),Z 
2- 

/ I + 7 \ 



-dE = x(z;Zo)F{E„,,Em) 



l+ZoJ 4TTDliz;zo)' 



where 



x(z;zo) : 



(A13) 



Em 



A-a„-T{E.z-A) 



E,„ Z)pp £(l+z)/(l+z„),z 



-dE 



Em 



E^-°'dE, 



(A14) 

is a function of the shape of the spectrum and the redshifts, 
similar to f{z;Zo)- Thus, a given heating rate defines an intrinsic 
luminosity limit for a source at a given redshift: 



imfcZo,^m) = 47r 



l+Z 
1+Zo 



a-2 



Dl(z;z„) 
x(z;zo) 



■qn 



(A15) 



With this, we may compute the heating rate as a function of 
qm, i.e., the heating rate associated with sources that produce a 
local heating larger than some limit: 



Q{z„\qm) = 



At:D\{z\Zo) 



dDp 

dz 

q4>B{z.,L)d\ogyQLdz 

log,|,L,„(c;;„4,„) 

2-a 



AnD\(z\Zo) 



dDp 



l+z 



dz \l+Zo 



(A16) 



X(z\Zo) 

4ttDI(z;zo) 



Ab [z;L„(z;zo,q„,),LM] dz. 



The heating from all gamma-ray blazars is obtained simply 
by setting q„, = 0. On the other hand, we may set implicitly 
via 

Qizo\qm)=QQ{zo;0), (A17) 

where Q ranges from to 1, yielding a heating rate limit at 
each redshift that we shall call qQ. From this, we may then 
obtain a number of contributing blazars: 



^fBjvizo; Q) ■■ 



, dDp 
dz 

X l>B [z\L„,{z;zo,qQ),J^M\ dz. (A18) 
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Figure 12. Radial distribution (x = D/Dt) of the sources that contribute to 
the total number (top) and heating rate (middle), and the limiting and effective 
luminosity (bottom) for a flux-limited sample of sources described by a fixed, 
cutoff power-law luminosity function in a static Euclidean Universe. Different 
colors coiTespond to different luminosity function power-laws (^ = —1, -0.8, 
-0.6, -0.4, —0.2, and shown in violet, blue, cyan, green, orange, and red), 
with the associated luminosity-weighted luminosity function shown explicitly 
in the inset. For reference, the radius at which D = Dpp is shown by the vertical 
dotted line. In comparison to the approximate L„, we used (black solid), we 
show the luminosity hmit when abso rption is included by the black dashed 
Hne, and the Leff defined in Equation )A25) ar e shown by the dot-dash lines. 
This may be compared directly with Figure ll4l 

This is shown for Q = 0.5 and Q = 0.75 in Figure [TT] by the 
black solid and long-dashed lines, respectively. While both are 
similar to the other measures of Mb below z ^ 1, above this 
redshift they fall much more rapidly. This is due to the shift 
of (j>B{z,L) towards higher luminosities, and thus the luminos- 
ity density (and therefore heating rate) becomes dominated by 
fewer, more luminous sources. Nevertheless, Mb is a rapidly 
increasing function of Q, as evidence by the fact that Mb.iv 
increases by approximately an order of magnitude when Q in- 
creases from 0.5 to 0.75. 

A.3. Location and Properties of the Sources Responsible for 
the Bulk of the Heating 

Armed with a definition of Mb, we can now address which 
class of high-energy gamma-ray blazars dominates the heat- 
ing rate. That is, we can assess whether the local heating is 
dominated by close, intrinsically dim objects or by distant, in- 
trinsically luminous sources. This is do ne sim ply by inspecting 
the integrands in Equations ( IA16I ) and ( IA18b . However, to in- 
terpret these, we will first build some intuition based upon an 
extremely simplified model, for which an analytical result is 
trivially obtained. 

A.3.1. Static Euclidean Universe 

We begin with a toy problem in which we consider a fixed 
power-law luminosity function in a static, Euclidean Universe. 
Specifically, we choose. 
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Figure 13. dloggfiA/logioL for a flux-limited sample of sources described 
by a fixed, cutoff power-law luminosity function in a static Euclidean Universe. 
Different colors correspond to different luminosity function power-laws, with 
the associated luminosity-weighted luminosity function shown explicitly in the 
inset. This may be compared directly with Figure [6] 



where 8(.ic) is the Heaviside function, is the overall normal- 
ization of the luminosity function, is a maximum luminos- 
ity (approximating a break), and ^ is an arbitrary constant. We 
will also assume that Dpp is independent of energy, i.e., there 
is some characteristic value for which t = D /D^^ and we may 
bring it out of the energy integral in the definition of q. As a 
consequence, a fixed limit in q„, corresponds directly to a fixed 
flux limit, F,„. 

In a Euclidean Universe the number of objects takes the par- 
ticularly simple form. 



\IE = I uiy^TTiy e ' " I (j){L) d\og^(^L 

/() "'logioi,,, 

V^m 



(A20) 



The flux limit, F„„ gives a luminosity limit of L,„ = 4ttD^F,„, 
where we have ignored the optical depth. Since we will be 
most concerned with how peaked the various integrands are at 
nearby distances, this is not a significant oversight (including it 
would serve to make them only more so). The definition of L„, 
implies a maximum distance as well, with = -^/L, /47rf„„ 
and therefore L,„/L^, = (D/D^,)- = x'. Thus, we have 



^fE = 



4-TT<f>^, 

In 10 JO 



l-(^ 



(A21) 
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From this we trivially obtain 

dAfs 1 dME 4TrDl 



dx^ 'x'' (l - 



dD Ds, dx 



In 10 
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(A22) 



providing some notion of the location of the most numerous 
sources. Typical values of ^ range from -1 to 0, in practice, 
and c/logA/fi/c/x is shown in the top panel of Figure [12] for a 
variety of choices of ^ within this range. 



The integral for Q is similarly simple, 

Qe= £/D47rDV«/«pp / mdlog^oL 
Jo "'iog,„L,„ 47rD^Dpp 



dx ^ -''-'pi 



lnlODpp7o 1+C 



and thus. 



dQE ^ i-x"^^^ ^_,/,^^ 

dD InlODpp 1-1-^ ^ 



(A23) 
(A24) 



giving an idea of the location of the sources responsible for the 
bulk of the heating. For a variety of ^, d\ogQE / dx is shown in 
the middle panel of Figure [12] 

Generally, we find that for all but the largest ^, dlogMs/dx 
and dlogQE /dx are peaked at small distances. In particu- 
lar, both are typically dominated by x < Xpp. The bottom 
panel of Figure [12] shows the flux-limited L,„ with (dashed) 
and without (solid) absorption included. Including absorption 
suppress the contributions at large x/xpp, forcing d\ogME / dx 
and dlogQE /dx to be even more strongly peaked at small dis- 
tances. 

Since L„, is a strong function of D, contributions from dif- 
ferent distances have different luminosity distributions. It is 
possible to roughly characterize this by defining a typical lu- 
minosity, Leff, associated with contributions at a given D: 



dQs/dD 
'dUE/dD 



e 1- 



25+2 



1+^ 1-x- 



2? 



(A25) 



With L„i, this is shown in the bottom panel of Figure [121 Gener- 
ally Leff is larger than L,„, and for small x substantially so. Thus, 
even for nearly flat 0(L) (^ ^ -1) the objects that contribute to 
the heating are not dominated by the numerous, intrinsically 
dim objects with luminosities L„,. 

Alternatively, we may perform the integral over D first, in 
which the flux limit implies a maximum distance to which a 
given object can be seen. Dm = y^L/4TrF„„ providing some in- 
sight into luminosity of the sources responsible for the heating. 
Doing so yields 



dQE 



t/logioL 




(A26) 



This is shown in Figure [13] for a variety of ^. In all cases 
the luminosity at the break in (piL), i.e., L*, contributes most 
significantly to the heating rate. However, the relative impor- 
tance of lower-luminosity objects does depend upon the lumi- 
nosity function; flat luminosity functions (i.e., ^ =-1) have 
many low-luminosity sources, and thus induce more broad 
d log QE/d logjo L. Nevertheless, it appears that the heating rate 
in our simple toy model is generally dominated by nearby ob- 
jects near the peak in the luminosity function. 

A.3.2. TeV Blazars in the Standard Cosmology 

We now return to the physically relevant case: heating due to 
TeV blazars with an evolving luminosity function in an evolv- 
ing Universe. Here we specify the distributions of the sources 
responsible for producing a fraction Q of the total heating rate. 
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Figure 14. Redshift distribution of the sources tiiat contribute to the total 
number (top-subpanel) and heating rate (middle-subpanel), and the limiting lu- 
minosity (bottom-subpanel) for a fractional Q-limited sample of objects. Dif- 
ferent colors coiTespond to different Zo, indicated by the left-most value of z 
for which each intersects the horizontal axis (z.o ranging from to 4), and are 
consistent with those used in Figure \6\ For reference, the redshifts at which 
f(z',Zii) = 1 and Lm are shown by the vertical and hoiizontal dotted lines, re- 
spectively. In comparison to L,„, we also show L^fi (as defined by Equation 
jA28)) by the dash-dot fines. Finally, tfie upper-rigfit panel sfiows (j)B(Zo,L) for 
eacfi of tfie redsfiifts for wfiicfi distributions are sfiown in tfie otfier panels, witfi 
corresponding colors (not e tfia t tfiis sfiows tfie same relative dynamic range as 
tfiat in tfie inset of Figures ll2l and may tfius be directly compared). 

In this case, we may immediately read off dAfs /dz and dQ/dz, 
the analog s of dNE/dD and dQE/dD, from Equations ( IA18b 
and ( IA16I ). respectively, yielding. 



dMs T dDp ~ r 1 

-{z;Zo) = 47rD;J {z;Zo) — <^B\z\L„{z\Zo-,qQ)-,I^M\ 
dz 
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dO J dDp 

-^(z-Zo) = ^^Dl(z-Zo)-^ 
dz dz 



\+z 

\+Zo 



2-Q 
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A^Di{z;zo) 

These are shown, normalized by their integrated values, for Zo 
ranging from 0-4 in Figure[T4]for Q = 0.5. In addition we show 
an analogously defined characteristic luminosity. 



^ 4ttDI{z;zo) dQ/dz 
~ X{z;Zo) dNs/dz 

1+z Afi [z;L,n{z;zo,qQ),LM\ 



(A28) 



1+Zo/ $B [z;L„,(z;Zo,?'Q),iM] 



The generic features of our static toy model are also appar- 
ent here. At all observer redshifts the heating rate is dominated 
by the nearest sources. At low Zo, where L(f)B is nearly flat, the 
number of objects is also heavily weighted towards nearby ob- 
jects, well within the redshift at which r = 1 . However, dNg /dz 
and dQ/dz evolve with observer redshift due to both, the intrin- 
sic evolution of (pB and the background Universe. As a conse- 
quence, by Zo ~ 0.5 the peak of dMB / dz has moved to the f = 1 



redshift, implying that z\ (where f{z\,Zo) = 1) is not a partic- 
ularly accurate estimate of the redshifts that contribute signifi- 
cantly to the heating. This is further supported by the high-Zo 
behavior of dJVs / dz, which once again is heavily weighted at 
redshifts inside of zi due to onset of the decline in the blazar 
population. 

The typical luminosities of objects responsible for the heat- 
ing also evolve. At Zo ^ these are roughly lO'^'^ergs"', rising 
to 3 X lO^'^ergs"' by Zo ^ 2. We also compute the heating rate 
per logarithmic decade in luminosity: 

dQ 



dlogi^L 



47tDI(z;z„) 



dDp 
~dz 



1+z 

XiZZo) 

4TrDl(z;zo) 



LMz.L), (A29) 



(where Zm is determined implicitly hy Lin{zm',Zo, 4 q) = ^-') shown 
in Figure |6] for a number of Zo- At Zo ~ the distribution of 
heating rates is a relatively broad function of L centered near 
10"*^ ergs. Until Zo ~ 2, as Zo grows dQ/dlog^^L becomes in- 
creasingly peaked and centered upon increasingly larger lumi- 
nosities. Above Zo 2 this trend reverses, though the distribu- 
tion of luminosities that contribute appreciably to the heating 
rate never becomes comparable to that in the present epoch. 
Thus, generally it appears that the heating is due predominantly 
to nearby objects with luminosities comparable to 10"*^ ergs"'. 
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